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Abstract 

In  this  dissertation,  we  investigate  the  minimum-time  orbital  transfer  problem  for 
spacecraft  with  steerable,  continuous  thrust  of  constant  magnitude.  The  optimal  control 
problem  is  developed  using  Euler- Lagrange  theory,  which  leads  to  the  optimal  control  law 
in  terms  of  the  Lagrange  multipliers  or  costates,  and  provides  the  differential  equations 
governing  these  costates.  Determination  of  the  costates  as  functions  of  time  through  nu¬ 
merical  solution  of  the  differential  equations  requires  initial  values  for  the  costates  which 
determine  the  initial  steering  angle  of  the  thruster.  It  is  well  known  that  finding  the  initial 
values  of  the  costates  is  the  most  difficult  part  of  solving  optimal  control  problems  of  this 
type.  The  standard  solution  technique  is  to  use  the  shooting  method  to  solve  a  boundary 
value  problem  in  which  the  initial  and  final  values  of  the  states  are  specified,  but  the  initial 
and  final  values  of  the  costates  are  unknown.  This  iterative  procedure  is  sensitively  depen¬ 
dent  on  the  initial  conditions  provided  for  the  costates.  This  research  has  developed  reliable 
approximate  models  for  the  initial  values  of  the  costates,  such  that  the  shooting  method 
will  always  converge  over  a  given  range  of  problem  parameters.  Employing  a  combination 
of  analytical  and  empirical  results,  the  optimal  initial  costates  are  modeled  as  functions 
of  the  problem  parameters  which  are  the  initial  thrust  acceleration,  A,  and  the  final  orbit 
radius,  i2,  in  canonical  units.  For  circle-to-circle,  coplanar  orbit  transfers,  these  approx¬ 
imate  initial  costate  models  lead  to  convergence  in  the  shooting  method  for  aU  practical 
values  of  A  and  i2.  In  addition,  the  models  lead  to  convergence  for  a  wide  range  of  other 
problems,  including  circle-to-hyperbola  transfers  and  non-coplanar  transfers.  To  counter 
the  extreme  sensitivity  to  small  changes  in  the  initial  costate  conditions,  a  dynamic  step 
limiter  is  introduced  which  improves  convergence  properties.  The  minimum-time  prob¬ 
lem  is  also  modeled  using  the  Kustaanheimo-Stiefel  (KS)  transformation,  and  the  optimal 
initial  costates  are  shown  for  comparison.  Several  numerical  examples  are  provided  for 
coplanar  and  non-coplanar  orbital  transfers  with  various  end  conditions. 
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L  Introduction 

1.1  Background 

The  fascinating  possibilities  of  space  flight  in  fact  and  fiction  have  inspired  many 
people  to  pursue  careers  in  the  astronautical  sciences.  As  a  result,  continual  advances 
are  being  made  in  the  science  of  spacecraft  design.  Propulsion  systems  in  particular  have 
improved  dramatically  since  the  first  black  powder  rockets  of  antiquity.  Robert  Goddard, 
also  inspired  by  the  potential  of  space  flight,  is  considered  the  father  of  modern  rocketry  for 
his  successful  experiments  using  liquid  fuels.  Today,  propulsion  research  continues  with  the 
development  of  high  efficiency  non-chemical  thrust  devices.  Thus,  there  are  three  common 
approaches  to  space  propulsion:  solid  fuel,  liquid  fuel,  and  non-chemical. 

Goddard  [13]  recognized  that  to  launch  his  rockets  in  the  most  fuel-efficient  manner 
possible,  he  would  need  to  solve  the  optimal  control  problem  using  the  calculus  of  varia¬ 
tions.  For  Goddard’s  experiments,  an  approximation  to  the  optimal  control  law  was  used 
instead.  Since  the  boost  phase  of  a  chemical  rocket  is  typically  of  short  duration,  an  ap¬ 
proximation  is  usually  considered  to  be  good  enough.  However,  a  full-scale  booster  rocket 
has  a  much  longer  thrusting  phase  than  Goddard’s  experiments.  Also,  the  thrust  duration 
of  a  non-chemical  propulsion  system  can  become  quite  lengthy.  In  either  case,  the  need 
for  optimal  control  solutions  becomes  much  greater  as  the  duration  of  the  thrusting  phase 
grows  large  compared  to  the  total  flight  time,  since  the  usual  approximations  become  less 
valid. 

Normally,  the  goal  of  the  optimization  problem  is  to  minimize  either  the  fuel  used  or 
the  time  taken  to  complete  a  mission.  The  problem  has  been  solved  [13]  for  both  quantities 
if  one  considers  the  velocity  changes  to  occur  instantaneously.  However,  there  is  stiU  no 
closed  form  solution  available  for  the  finite  duration  thrust  case.  Numerical  methods  are 
typically  used  to  determine  an  optimal  thrust  program  to  meet  boundary  conditions  for 
position  and  velocity. 
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Non-chemical  thrust  devices  are  appealing  because  they  typically  have  high  propul¬ 
sive  efficiencies,  measured  by  the  specific  impulse.  A  highly  efficient  propulsive  system  is 
very  desirable  for  an  orbiting  spacecraft,  simply  because  there  is  rarely,  if  ever,  an  opportu¬ 
nity  to  refuel  on  orbit.  Although  finite  thrust  devices  such  as  electrical  arcjets  use  less  fuel 
than  chemical  thrusters  for  a  given  change  in  spacecraft  velocity,  they  normally  produce 
very  small  thrust  levels.  For  this  reason,  finite  thrust  propulsion  systems  need  to  operate 
continuously  for  extended  periods  to  accomplish  orbital  maneuvers. 

L2  Problem  Statement 

When  a  spacecraft  is  being  accelerated  by  a  thruster  for  significant  portions  of  the 
planned  trajectory,  the  effect  can  not  be  considered  instantaneous,  or  “impulsive.”  Thus, 
the  orbital  path  will  not  be  Keplerian.  The  optimal  magnitude  and  direction  of  the  thrust 
must  then  be  found  as  a  function  of  time  to  meet  mission  objectives.  It  is  possible  to 
find  this  function,  the  optimal  control  law,  by  using  Euler- Lagrange  theory,  which  will  be 
discussed  in  Chapter  4.  However,  knowledge  of  the  optimal  control  law  is  not  sufficient  to 
solve  the  problem  of  meeting  desired  end  conditions,  because  Euler- Lagrange  theory  intro¬ 
duces  adjoint  variables  which  must  be  initialized.  These  variables,  also  known  as  Lagrange 
multipliers  or  costates,  are  difficult  to  initialize  because  there  is  insufficient  information 
from  the  boundary  conditions  to  specify  their  initial  values.  Without  this  information,  it  is 
not  possible  to  propagate  the  differential  equations  that  govern  the  behavior  of  the  states 
and  costates.  Typically,  the  initial  values  of  the  costates  are  guessed,  then  an  attempt  is 
made  to  solve  the  boundary  value  problem  by  refining  the  guesses  in  some  automated  way. 
Due  to  the  sensitivity  of  the  costates  to  errors  in  the  initial  conditions,  poor  guesses  may 
preclude  any  hope  of  convergence  to  the  desired  end  conditions.  Thus,  the  problem  is  to 
find  the  initial  values  of  the  costates  that  will  lead  to  the  desired  final  orbital  conditions, 
in  the  minimum  time. 

L3  Research  Goal 

The  goal  of  this  research  is  to  provide  insight  into  the  selection  of  initial  values  for 
the  Lagrange  multipliers,  leading  to  reasonable  certainty  of  convergence  for  the  boundary 
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value  problem  of  a  spacecraft  under  continuous  thrust  with  fixed  end  conditions.  To  gain 
this  insight,  both  analytical  and  empirical  means  will  be  used  to  model  the  optimal  initial 
costates  as  functions  of  the  problem  parameters. 

L4  Thesis  Outline 

Chapter  2  provides  a  literature  review  related  to  the  problem  of  optimal  control 
for  impulsive  and  continuous-thrust  orbital  maneuvers.  In  Chapter  3,  the  equations  of 
motion  for  a  spacecraft  influenced  by  gravity  and  continuous  thrust  are  derived  in  several 
coordinate  systems.  Chapter  4  starts  with  a  presentation  of  Euler- Lagrange  theory,  which 
is  then  used  to  develop  the  optimal  control  law  and  costate  equations  in  each  of  the 
coordinate  systems.  The  shooting  method  is  also  discussed  in  Chapter  4.  Chapter  5 
contains  an  analysis  of  the  optimal  initial  costate  locus,  which  is  used  to  initialize  the 
shooting  method.  Numerical  examples  of  optimal  continuous-thrust  transfers  using  the 
costate  locus  analysis  are  presented  in  Chapter  6.  Finally,  a  summary  of  this  research  is 
presented  and  conclusions  are  discussed  in  Chapter  7. 
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IL  Literature  Review 


There  are  hundreds,  if  not  thousands,  of  papers  to  be  found  in  the  literature  on 
optimal  space  maneuvers.  For  example.  Bell  [5]  cites  160  articles  in  his  survey  of  pub¬ 
lished  work  on  optimal  space  trajectories.  The  two  most  commonly  addressed  issues  are 
minimum-fuel  transfers  and  minimum-time  transfers.  For  a  spacecraft  under  continuous 
thrust  with  no  throttling,  the  resulting  solutions  wiU  be  the  same.  Without  throttling,  the 
mass  flow  rate  of  propellant  is  constant,  so  if  time  is  minimized  then  the  fuel  consumed 
will  be  minimized  as  well. 

Depending  on  the  design  of  the  propulsion  system,  the  velocity  change  may  occur  in 
a  very  short  time,  or  in  a  very  long  time  compared  to  the  period  of  the  desired  final  orbit. 
Short  thrusts  are  treated  as  impulsive^  and  long  thrusts  are  modeled  as  continuous  effects 
for  finite  durations. 

2 A  Optimal  Impulsive  Maneuvers 

One  of  the  earliest  definitive  works  on  optimized  impulsive  maneuvers  is  by  Law- 
den  [13].  He  posed  the  minimum-fuel  space  trajectory  as  a  Mayer  problem  [13],  and 
sought  solutions  using  variational  calculus  methods  and  Lagrange  multipliers.  Lawden 
treated  the  Lagrange  multipliers  as  components  of  a  vector,  which  he  called  the  “primer 
vector.”  The  behavior  of  the  primer  vector  gave  the  optimal  directions  for  impulsive  ma¬ 
neuvers,  and  thus  solved  the  optimization  problem  for  impulsive  thrust.  This  work  also 
verified  Hohmann’s  result  for  a  minimum-fuel  impulsive  orbital  transfer.  Lawden’s  book 
is  commonly  referenced  in  contemporary  literature  and  serves  as  a  starting  point  for  much 
of  the  work  that  follows. 

2,2  Optimal  Continuous  Thrust  Maneuvers 

Although  the  optimization  of  impulsive  transfers  yielded  a  direct  solution  [13],  none 
has  been  found  for  the  continuous-thrust  case.  This  problem  may  be  solved  numeri¬ 
cally,  and  many  examples  of  this  are  to  be  found  in  the  literature  [5].  Optimization  of 
a  continuous-thrust  trajectory  involves  the  simultaneous  solution  of  an  optimal  control 
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problem  and  a  boundary  value  problem.  The  initial  and  final  states  are  normally  given, 
but  there  is  usually  no  information  available  for  the  initial  values  of  the  Lagrange  multipli¬ 
ers.  This  presents  quite  a  problem,  since  the  optimal  control  law  is  often  a  function  of  the 
Lagrange  multipliers  which  must  be  initialized  for  numerical  integration.  The  usual  ap¬ 
proach  is  to  make  an  educated  guess  for  the  initial  values,  then  update  them  by  solving  the 
boundary  value  problem.  Trussing  [17],  Broucke  [6]  and  others  have  recast  the  boundary 
value  problem  in  terms  of  other  variables,  but  the  initial  values  of  these  must  be  guessed 
and  refined  as  well.  Trussing  [17]  incorporates  the  second  derivative  of  the  primer  vector 
into  a  fourth  order  dynamics  equation,  thus  eliminating  the  control  variables.  Once  this 
is  accomplished,  four  constants  of  integration  must  be  iterated  to  find  the  correct  optimal 
trajectory.  Broucke  [6]  expresses  the  Lagrange  multipliers  as  functions  of  new  auxiliary 
variables,  and  graphically  examines  the  behavior  of  the  new  variables.  Tines  [15]  and  Red¬ 
ding  and  Breakwell  [18]  have  suggested  using  the  results  of  optimal  impulsive  maneuvers 
to  serve  as  an  initial  guess  for  the  continuous-thrust  case.  However,  this  method  produces 
poor  results  for  small  values  of  continuous  thrust,  particularly  if  there  are  no  coasting  arcs 
used.  Thus,  there  are  no  models  or  techniques  in  the  literature  to  provide  good  estimates 
for  the  initial  Lagrange  multiplier  values  for  the  continuous-thrust,  minimum-time  orbit 
transfer  problem. 

Closed  form  non-optimal  solutions  have  been  found  for  spacecraft  trajectories  where 
special  assumptions  are  made  about  the  control  law.  If  the  thrust  vector  is  directed  either 
radially  from  the  attracting  center  or  tangentially  to  the  orbital  path,  it  is  possible  to  inte¬ 
grate  the  equations  of  motion  analytically.  Battin  [4]  (section  8.8)  presents  results  for  the 
time  to  reach  escape  velocity  and  the  number  of  revolutions  for  both  thrust  assumptions. 

Assumptions  about  the  thrust  magnitude  will  also  allow  closed  form  non-optimal 
solutions  through  the  method  of  averaging  [1,  26].  If  the  thrust  level  is  small  enough,  there 
is  only  a  small  change  in  semi-major  axis  or  eccentricity  for  a  single  orbital  revolution. 
Then,  a  correction  is  made  to  the  semi-major  axis  at  the  completion  of  each  revolution. 
These  approximations  are  reasonable  for  orbital  transfers  that  require  roughly  ten  or  more 
revolutions  to  complete  [1,  26].  Using  these  assumptions,  it  is  possible  to  solve  for  the 
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trajectory  analytically.  The  thrust  is  directed  tangentially,  which  is  also  perpendicular  to 
the  orbit  radius  since  the  eccentricity  is  assumed  to  be  zero  for  individual  revolutions. 

In  another  approach  [11],  finite  difference  equations  are  used  in  place  of  the  exact 
differential  equations  of  motion.  Then,  a  choice  must  be  made  for  a  step  or  mesh  size  in  the 
search  space.  By  refining  the  mesh,  the  solution  may  approach  the  optimal  trajectory.  A 
method  known  as  ‘‘differential  inclusion”  [19,  8]  also  uses  finite  difference  equations.  These 
methods  can  be  very  efficient.  However,  it  is  difficult  to  guarantee  that  the  converged 
solution  is  the  desired  optimal  path  since  the  differential  equations  governing  the  costates 
are  not  used.  These  methods  have  gained  in  popularity  because  of  their  inherent  robustness 
and  the  increasing  power  of  digital  computers. 

Another  numerical  method  that  has  been  used  with  success  for  the  minimum-fuel 
problem  is  hybrid  non-linear  programming,  or  HNLP  [27].  In  this  case,  the  cost  function 
is  evaluated  directly  while  the  transversality  conditions  are  satisfied  implicitly.  HNLP 
combines  the  advantages  of  using  costate  equations  and  the  simplicity  of  directly  evaluating 
the  cost  function.  This  method  can  be  made  more  robust  than  propagating  the  exact  state 
and  costate  equations,  but  the  performance  depends  on  an  optimal  choice  of  additional 
variables  and  constraints. 

For  circle-to-circle  coplanar  orbital  transfers,  the  minimum  time  of  flight  may  be 
derived  from  the  accumulated  velocity  change  on  the  trajectory  [2].  It  is  possible  to  display 
the  optimal  accumulated  velocity  change  in  graphical  form  as  a  function  of  constant  thrust 
level,  ratio  of  final  to  initial  orbit  radius,  and  mass  propellant  fraction  [2].  In  this  way,  a 
wide  range  of  possible  cases  may  be  represented  through  the  use  of  universal  variables.  To 
produce  the  graphical  results,  many  different  cases  must  be  solved  numerically  to  allow  for 
interpolation.  Although  linear  interpolation  from  a  graph  wiU  not  provide  great  precision, 
it  does  show  general  trends  for  mission  design.  In  particular,  a  graph  of  the  number  of 
revolutions  for  the  optimal  path  versus  the  logarithm  of  thrust  magnitude  shows  a  distinct 
change  in  trajectory  characteristics  at  integer  values  of  revolutions  [2], 


2-3 


2,3  Summary 

Many  different  approaches  have  been  used  to  solve  the  optimal  continuous-thrust 
orbit  transfer  problem.  While  methods  that  propagate  the  exact  equations  of  the  states 
and  costates  will  guarantee  optimality,  they  are  not  robust  due  to  sensitivity  to  initial 
conditions.  Direct  methods  that  simply  evaluate  the  cost  function  with  approximate  finite 
difference  equations  are  typically  robust,  but  are  not  optimal.  Finally,  exact  solutions  to 
approximations  of  the  equations  of  motion  or  control  law  also  sacrifice  optimality.  An  ideal 
solution  technique  would  have  both  guaranteed  optimality  and  robustness,  but  none  exists 
in  the  literature  for  a  wide  range  of  spacecraft  design  parameters  and  orbital  boundary 
conditions.  In  particular,  there  are  no  models  or  techniques  in  the  literature  to  provide  good 
estimates  for  the  initial  Lagrange  multiplier  values  for  the  continuous-thrust,  minimum¬ 
time  orbit  transfer  problem,  based  on  problem  parameters. 
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III.  Continuous- Thrust  Spacecraft  Dynamics 

In  this  chapter,  the  equations  of  motion  for  a  spacecraft  under  the  influence  of  grav¬ 
ity  and  continuous  thrust  are  derived  for  a  variety  of  coordinate  systems.  The  resulting 
differential  equations  will  be  used  for  the  optimal  control  formulations  of  Chapter  4.  The 
acceleration  due  to  constant,  non-throttleable  thrust  is  A,  the  gravitational  parameter  is 
fjL.  The  length  of  the  position  vector  is  r,  and  the  final  desired  value  of  r  will  be  given  by  R. 
The  spacecraft’s  initial  mass  is  mo,  and  the  mass  flow  rate  is  m.  The  thrust  acceleration 
A  may  or  may  not  be  a  function  of  time,  depending  on  the  value  of  m,  as  will  be  shown  in 
the  next  section.  In  canonical  units  [3],  the  gravitational  constant  /x  is  unity  regardless  of 
the  system  under  consideration  as  long  as  the  initial  radius  is  defined  to  be  one  distance 
unit,  (DU),  and  the  initial  circular  velocity  at  that  radius  is  one  distance  unit  per  time 
unit,  (DU/TU).  Also,  the  initial  spacecraft  mass  is  one  mass  unit  (MU). 

3. 1  Equations  of  Motion  in  Three  Dimensions 


*  y’ 


X 

Figure  3.1  Problem  Geometry  in  Three  Dimensions 

The  equations  of  motion  are  most  easily  expressed  using  an  inert iaUy  fixed,  right- 
handed  Cartesian  system  o{  and  coordinates  as  shown  in  Figure  3.1.  The  and  y' 
axes  shown  in  Figure  3.1  are  parallel  to  the  x  and  y  axes,  respectively.  The  initial  circular 
orbit  lies  in  the  x,  y  plane,  and  the  initial  position  is  at  x  =  1,  y  =  0,  2:  =  0.  The  thrust 
direction  must  be  defined  with  (at  least)  two  angles.  The  angle  a  lies  in  the  x,  y  plane, 
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and  is  measured  from  the  x  axis  in  the  positive  direction.  The  angle  j3  is  measured  “up” 
from  the  y  plane,  and  lies  in  the  plane  containing  the  thrust  vector  and  the  z  axis.  The 
magnitude  of  the  position  vector  is  r  =  y/x’^  + 

The  equation  of  motion  for  a  rocket  [25]  is  as  follows: 

T>Fext  +  =  mf  (3.1) 

where  '^F^xt  is  the  sum  of  the  external  forces  acting  on  the  spacecraft,  is  the  exhaust 
velocity  of  the  propellant,  f  is  the  position  vector,  and  f  is  the  spacecraft  total  acceleration. 
The  second  term  of  Equation  (3.1)  represents  the  thrust,  and  the  magnitude  of  the  thrust 
acceleration  is  A  =  T/{mo+mt)  where  T  is  the  constant  thrust  magnitude  of  the  propulsion 
system,  and  t  is  the  time.  If  m  =  0,  then  A  is  equal  to  a  constant.  The  only  external 
force  we  will  consider  is  that  of  gravity  from  a  single  point  source,  and  Newton’s  Law  of 
Gravitation  may  be  used  to  express  this  as  follows: 

f,  =  (3-2) 

The  two-body  assumption  is  used  to  simplify  the  equations  of  motion  enough  to  allow 
closed  form  solutions  which  are  presented  in  Chapter  5.  Other  gravitational  models  may  be 
considered,  but  the  majority  of  missions  of  interest  are  primarily  influenced  by  two-body 
effects  and  thrust.  Therefore,  any  results  obtained  with  the  two-body  assumption  may 
be  used  as  a  starting  point  for  a  great  number  of  more  complex  orbital  problems.  It  may 
even  be  possible  to  approach  the  restricted  three-body  problem  by  treating  the  mass  of  the 
third  body  as  a  parameter,  and  using  two-body  results  as  initial  estimates  for  small  values 
of  the  mass  parameter.  However,  that  subject  is  beyond  the  scope  of  this  research.  The 
equations  of  motion  are  as  follows,  where  the  thrust  terms  are  the  components  of  mV^: 


X  = 

—{lifr^)x  H-  A  cos  /3  cos  a 

(3.3) 

y  = 

—{fi/r^)y  +  A  cos  f3  sin  a 

(3.4) 

z  = 

—{ii/t^)z  +  a  sin  (3 

(3.5) 
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3.2  Equations  of  Motion  in  Two  Dimensions 


For  the  case  where  the  motion  remains  in  the  initial  orbital  plane,  the  geometry 
simplifies  to  that  of  Figure  3.2.  Since  /3  =  0  =  Equations  (3.4)-(3.5)  reduce  to: 


— (/x/r^)x  +  A  cos  a 

(3.6) 

— (/x/r^)y  +  A  sin  a 

(3.7) 

The  two-dimensional  equations  of  motion  may  also  be  expressed  in  polar  coordinates 
which  include  r  as  the  scalar  distance  from  the  attracting  center,  u  as  the  time  rate 
of  change  of  r,  and  v  as  the  velocity  component  perpendicular  to  u  directed  along  the 
spacecraft  horizon.  The  polar  thrust  angle,  cj)^  is  measured  clockwise  (“up”)  from  the 
spacecraft  local  horizontal,  as  shown  in  Figure  3.2.  This  results  in  the  following  differential 
equations: 


T  =  u  (3.8) 

2 

ii  =  - - +  Asincf)  (3.9) 

T 
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V 


(3.10) 


uv  ^ 

— - \-  A  cos 

T 
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3,8  Equations  of  Motion  under  the  KS  Transformation 

The  Kustaanheimo-Stiefel  (KS)  transformation  is  intended  to  regularize  the  equa¬ 
tions  of  motion  in  the  problem  of  two  bodies  [21].  The  purpose  of  the  regularization  is  to 
reduce  numerical  integration  dilRculties  when  r  is  small,  by  placing  the  inverse  of  r  into 
a  term  that  represents  the  constant  angular  momentum  magnitude  of  a  two-body  orbit. 
This  term  premultiplies  the  state  variables  Ui  and  t^2  in  Equations  (3.15)  and  (3.16),  but  it 
will  not  remain  constant  with  the  influence  of  thrust.  When  this  transformation  is  used  in 
conjunction  with  a  change  of  independent  variable,  the  equation  of  motion  in  two  dimen¬ 
sions  has  the  form  of  a  harmonic  oscillator  [20].  This  allows  for  simple  analytical  solutions, 
which  may  be  perturbed  by  other  forces  such  as  a  third  body  or  a  propulsion  system.  In 
Chapter  4,  we  wiU  apply  Euler- Lagrange  theory  to  the  two-dimensional  regularized  equa¬ 
tions  of  motion  to  solve  the  minimum-time  problem  for  a  coplanar  orbital  transfer  under 
continuous  thrust.  The  purpose  of  this  development  is  to  compare  the  costate  behavior 
with  the  Cartesian  case  in  Chapters  5  and  6. 

The  equations  of  motion  for  a  two-dimensional  orbit  are  as  follows: 

(3.11) 

(3.12) 


Using  the  KS  transformation  for  two  dimensions,  the  coordinates  x  and  y  are  replaced 
by  Ui  and  U2  through  the  following  relationship: 

(til  +  iu2f  =  X  +  iy  (3.13) 


The  independent  variable  t  is  replaced  by  the  fictitious  time  s  with  the  following  differential 
equation: 


dt 

ds 


=  r 


(3.14) 
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These  transformations  lead  to  the  regularized  equations  of  motion  [20]: 


Uo  = 


(3.15) 

(3.16) 


The  primes  indicate  differentiation  with  respect  to  s,  stnd  r  =  uj  +  uj.  The 

symbol  indicates  an  inner  product.  A  thrust  model  may  be  added  as  follows: 


Ui  +  ^Ar^/^cos7 


U2  +  -Ar^/^sin7 


(3.17) 


(3.18) 


where  A  is  the  magnitude  of  the  thrust  acceleration.  This  is  an  original  thrust  model  which 
is  consistent  with  the  transformation  given  in  Reference  [20],  Equation  (9.26),  but  here  the 
thrust  angle  7  has  been  defined  in  the  u  coordinate  system  for  simplification  with  no  loss 
of  generality.  The  relationship  between  the  inertial  Cartesian  thrust  angle,  a,  and  the  KS 
thrust  angle,  7,  is  as  follows: 


cos 7  =  T  COSO!  +  U2  sin  a)  (3.19) 

sin  7  =  sin  o  —  ^2  cos  o)  (3.20) 


At  the  initial  time,  Ui  =  1  and  U2  =  0,  and  r  =  1.  Therefore,  7(0)  =  o(0).  Defining  Vi  = 
the  equations  of  motion  and  differential  constraints  may  be  expressed  as  five  first-order 
differential  equations: 


T 

(3.21) 

Vl 

(3.22) 

V2 

(3.23) 
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(3.24) 


V,  = 


Ui  +  cos  7 


+  si 


sm7 


(3.25) 


S,4  Summary 

The  equations  of  motion  of  a  spacecraft  under  the  influence  of  gravity  and  continu¬ 
ous  thrust  have  been  presented  in  four  different  coordinate  frames.  In  three  dimensions, 
inertial  Cartesian  coordinates  are  used.  In  two  dimensions,  inertial  Cartesian  or  polar 
coordinates  are  used.  Finally,  the  equations  of  motion  are  modified  by  the  KS  transfor¬ 
mation.  An  original  thrust  model  is  presented  that  is  simplified  compared  to  the  model 
given  in  Reference  [20].  These  sets  of  equations  will  be  used  in  Chapter  4  to  derive  optimal 
control  formulations  for  the  minimum-time  orbital  transfer  problem. 
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IV.  Optimal  Control  Formulations  and  Solutions 

This  chapter  begins  with  a  standard  development  of  the  variational  calculus  approach 
to  the  minimum-time  trajectory  optimization  problem  [7].  Then,  the  resulting  relationships 
are  applied  to  the  equations  of  motion  for  a  spacecraft  under  continuous  thrust  in  several 
coordinate  systems,  as  derived  in  Chapter  3.  Finally,  the  shooting  method  is  presented 
as  a  means  to  solve  the  two-point  boundary  value  problem  of  transferring  from  one  orbit 
to  another,  using  the  optimal  differential  equations  for  the  states  and  costates.  There  are 
two  original  presentations  in  this  chapter:  the  application  of  Euler- Lagrange  theory  to 
the  equations  of  motion  under  the  KS  transformation,  and  a  dynamic  step  limiter  which 
improves  the  convergence  properties  of  the  shooting  method. 

Optimization  problems  come  in  many  forms,  but  the  usual  goal  is  to  minimize  (or 
maximize)  the  value  of  some  desired  quantity.  Sometimes,  the  answer  we  seek  is  a  function 
of  time  that  will  minimize  this  quantity.  Suppose  it  is  desired  to  get  from  one  position  to 
another  in  a  minimum  time.  The  solution  to  this  problem  is  a  path,  or  set  of  directions 
as  a  function  of  time.  For  example,  a  map  with  highlighted  roads  from  an  automobile 
club  could  be  considered  an  optimal  path  to  minimize  travel  time  between  cities.  A  better 
example  might  be  a  program  for  a  road  raUy,  in  which  the  drivers  are  expected  to  arrive  at 
intermediate  checkpoints  at  specific  times.  In  this  case,  the  goal  would  be  to  maximize  a 
score,  rather  than  to  minimize  the  time.  Either  way,  the  set  of  instructions  specifies  quan¬ 
tities  like  position,  speed,  and  direction.  To  follow  an  optimal  path  as  in  these  examples,  a 
driver  must  control  the  vehicle  by  steering,  accelerating,  braking,  and  so  on.  These  control 
actions  must  be  done  in  exactly  the  right  order  at  the  right  time  for  success.  Thus,  a  set 
of  control  instructions  or  control  law  is  the  solution  that  we  ultimately  require. 

Euler- Lagrange  theory  provides  an  analytical  method  to  find  the  control  law  for  a 
variety  of  path  optimization  problems  through  the  calculus  of  variations  [7].  To  implement 
this  theory,  we  need  to  define  a  quantity  to  be  minimized  (or  maximized),  and  to  determine 
any  constraints  that  affect  the  optimal  path.  For  a  problem  in  spaceflight,  the  vehicle 
must  move  according  to  Newton’s  laws  of  motion  under  the  influence  of  gravity  and  the 
propulsion  system.  We  will  not  consider  other  perturbations  or  relativistic  effects  on  the 
spacecraft.  The  equations  of  motion  may  be  treated  as  constraints,  since  the  vehicle  has 
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no  other  choice  than  to  obey  them.  These  constraints  may  be  written  as  a  set  of  first  order 
dilferential  equations: 

x{t)  ^  (4.1) 

If  we  wish  to  minimize  or  maximize  a  quantity  a  cost  functional  is  defined  as: 


where  tj  is  the  final  time.  The  integral  form  is  chosen  so  that  the  integrand  may  be 
identified  as  the  Lagrangian  (£),  which  wiU  be  used  in  the  calculus  of  variations  approach. 
The  constraints  may  be  added  to  the  Lagrangian  without  changing  the  value  of  the  cost 
function,  if  they  are  appended  in  the  following  way: 

J  ^  (c  +  x'^  (/-  f))  dt  (4.3) 

Notice  that  the  quantity  (^f  —  is  identically  equal  to  zero  if  the  state  vector,  Xj  solves  the 
differential  equations,  /,  by  the  definition  of  the  equations  of  motion  given  above.  Thus, 
we  have  added  nothing  to  the  cost  function,  as  long  as  the  constraints  are  obeyed.  The 
variable  A  is  a  vector  of  scalar  valued  functions  known  as  Lagrange  multipliers.  Strictly 
speaking,  the  vector  is  a  representation  of  linear  functionals  [14].  The  Lagrange  multipliers 
are  also  known  as  costates,  since  there  will  be  one  element  of  the  Lagrange  multiplier  vector 
associated  with  each  element  of  the  state  vector.  The  behavior  of  the  costates  turns  out 
to  be  very  important  for  defining  the  optimal  control  law,  as  will  be  shown. 

The  cost  function  is  now  dependent  on  the  states,  £,  and  the  controls,  fT,  which  are 
contained  in  /.  An  infinitesimal  variation  away  from  the  minimum  point  wiU  not  increase 
the  value  of  the  cost,  since  the  slope  is  zero  there.  The  variation  operation  is  similar  to 
taking  a  derivative,  but  it  provides  information  about  the  relationships  of  many  variables 
at  the  same  time.  First,  define  a  scalar  valued  function  the  Hamiltonian: 

=  £  +  F/  (4.4) 
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The  first,  total  variation  of  the  cost  function  for  either  fixed  or  free  final  conditions 
and  fixed  initial  conditions  is  given  by  [12]: 


6x  + 


dt 


(4.5) 


We  would  like  to  be  able  to  make  small  variations  in  the  control  and  final  time  without 
changing  the  value  of  the  cost  function.  For  this  to  happen,  the  variation  of  J  must  be 
zero  for  arbitrary  Stj  and  Su.  The  states  must  be  free  to  vary  except  at  the  endpoints, 
so  6x  is  also  arbitrary  in  the  integrand.  If  some  of  our  final  conditions  are  not  specified, 
such  as  the  final  radius,  then  the  first  term  may  have  non-zero  elements  for  6x{tj).  Since 
we  assume  the  variation  of  the  cost  function  to  be  identically  zero,  the  final  Lagrange 
multiplier  vector  would  have  as  many  zero  components  as  necessary  for  the  first  term  to 
vanish.  As  mentioned  above,  this  reasoning  will  also  apply  if  the  initial  states  are  free  and 
the  final  states  are  specified  [12,  7]. 

The  second  term  in  the  variation  requires  that  the  final  value  of  the  Hamiltonian, 
7f(ty),  is  equal  to  zero  if  6tf  is  free  to  vary.  This  is  known  as  the  transversality  condition  [7]. 
For  minimum-time  problems,  the  final  time  must  be  free  to  vary,  so  the  final  value  of  the 
Hamiltonian  must  be  zero.  This  condition  may  be  met  through  proper  scaling  of  the 
Lagrange  multipliers.  Scaling  may  also  be  used  to  make  the  first  term  in  the  variation  zero 
as  weU.  The  equations  of  motion  for  the  states  and  costates  in  our  problem  will  not  be 
affected  by  scaling  the  Lagrange  multipliers,  which  gives  us  some  freedom  to  choose  the 
scale  factor.  This  is  because  the  Lagrange  multipliers  appear  in  the  equations  of  motion 
for  the  state  variables  as  ratios  of  each  other,  as  is  shown  in  the  following  sections. 

The  integral  term  in  the  variation  must  also  be  zero  for  arbitrary  values  of  6x  and 
6u.  Thus,  for  the  integrand  to  be  zero,  we  have  the  Hamilton- Jacobi  equations  [11]: 


^  ^  nT 

du 

(4.6) 

dn  f 

dx 

(4.7) 
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Equation  (4.6)  is  the  optimality  condition,  which  states  that  the  variation  of  the  Hamil¬ 
tonian  with  respect  to  the  control  should  be  zero  on  the  optimal  path.  Equation  (4.7) 
provides  a  set  of  first  order  differential  equations  that  govern  the  behavior  of  the  Lagrange 
multipliers.  These  are  the  costate  equations,  which  may  be  integrated  along  with  the  state 
equations  through  the  time  interval. 

In  the  calculus  of  variations,  the  Legendre-Clebsch  condition  is  used  to  check  if  the 
control  u  minimizes  or  maximizes  the  Hamiltonian.  For  minimization,  we  have  [7]: 


which  states  that  is  a  positive  semi-definite  matrix  throughout  the  time  interval. 

For  maximization,  we  have: 

S<0  (4.9) 


which  states  that  is  a  negative  semi-definite  matrix  throughout  the  time  interval. 

Thus,  we  may  check  for  either  occurrence.  If  the  second  partial  derivative  of  the  Hamilto¬ 
nian  with  respect  to  the  control  is  equal  to  zero  or  a  zero  matrix,  then  we  have  a  singular 
arc  [12],  and  further  investigation  is  necessary  to  determine  the  nature  of  the  critical  path. 
However,  the  partial  derivatives  of  the  Hamiltonian  are  functions  of  time,  and  numerical 
experience  has  shown  that  the  second  partial  derivative  is  extremely  unlikely  to  remain 
exactly  equal  to  zero  throughout  a  trajectory.  This  situation  would  amount  to  a  case  in 
which  the  optimal  solution  does  not  depend  on  the  constraint  equations  in  the  Hamilto¬ 
nian,  as  may  be  seen  from  control  laws  which  are  developed  in  the  following  sections.  The 
Legendre-Clebsch  condition  provides  the  necessary  conditions  for  a  minimizing  path  [7]; 
however,  to  have  both  necessary  and  sufficient  conditions  for  a  minimum,  one  must  make 
certain  that  there  are  no  singular  arcs  encountered  along  the  path.  In  the  problem  we  are 
investigating,  the  second  partial  derivative  is  positive  or  negative  semi-definite  in  all  cases 
of  interest. 


The  optimality  condition  and  the  costate  equations  are  a  very  powerful  result  of 
Euler-Lagrange  theory.  This  result  may  be  used  to  determine  the  optimal  control  law  for 


4-4 


a  spacecraft  under  continuous  thrust.  However,  there  is  a  significant  difficulty  inherent  in 
this  formulation  which  is  addressed  in  this  research. 

Although  the  costate  equations  may  be  derived  from  the  above  relationships,  they 
must  be  initialized  to  begin  a  numerical  integration  procedure.  We  may  choose  the  initial 
and  final  conditions  for  the  physical  states,  but  there  is  no  guaranteed  way  to  determine  the 
correct,  optimal  boundary  conditions  for  the  costates.  In  addition,  the  costate  equations 
tend  to  be  very  sensitive  to  initial  conditions  in  practice.  They  are  also  just  as  sensitive 
when  choosing  final  conditions  for  backwards  integration,  as  in  a  forward/backward  sweep 
approach  [7].  The  shooting  method,  which  is  described  in  a  later  section,  depends  on 
“reasonable”  choices  of  initial  Lagrange  multiplier  values.  If  they  are  too  far  away  from 
the  correct  values,  the  shooting  method  will  fail. 


4.1  Optimal  Control  in  Three  Dimensions 

Using  the  three-dimensional  equations  of  motion  derived  in  Chapter  3,  Equations 
(3.3)-(3.5),  we  may  form  the  following  variational  Hamiltonian: 

1-1  =  1  +  Xx^  +  XyP  -h  XzZ  +  Aj;  [—{fji/r^)x  +  A  cos  /?  cos  a] 

+  Xy  [— (/i/r^)x  +  A  cos /3  sin  a]  +  Ai  [—(p,/r^)x  -b  A  sin/?]  (4-10) 


The  optimality  condition  leads  to  two  control  laws  for  the  thrust  vector  angles  a  and  /?: 


tan  a 

tan  /? 


Xx 

-A, 

—  (Ai  cos  a  +  Xy  sin  a) 


(4.11) 

(4.12) 


This  choice  of  sign  will  guarantee  the  necessary  conditions  for  a  minimum  with  respect  to 
each  of  the  control  angles,  as  is  shown  next  using  the  Legendre- Clebsch  condition.  Since 
there  are  two  control  angles,  the  second  partial  derivative  of  the  Hamiltonian  will  be  a 
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2x2  matrix: 


d^H/da^  d'^'H.lid^da) 

dmi{dadi3)  d^n/d/s^ 

Starting  with  the  first  diagonal  term,  we  have: 


d^n 

du 


d^n 

da^ 


—XiA  cos  (3  cos  a  —  cos  /3  sin  a 


(4.13) 


(4.14) 


Using  Equations  (4.11)  and  (4.12)  to  eliminate  a  and  (3  from  the  right-hand  side,  it  is 
possible  to  determine  the  sign  of  the  second  partial  derivative  of  the  Hamiltonian  with 
respect  to  the  control  angle  a: 


^A|  +  A?  +  A| 


>  0 


(4.15) 


The  second  diagonal  term  is  given  by  the  second  partial  derivative  of  the  Hamiltonian  with 
respect  to  the  control  angle 


d^n 

5/32 


Aii4  cos  a(—  cos/3)  +  XyA  sin  a{—  cos /3)  —  XiA  sin  /3 


(4.16) 


As  before,  Equations  (4.11)  and  (4.12)  may  be  used  to  determine  the  sign  of  the  second 
partial  derivative  with  respect  to  /3: 


d^n  .  a|  +  a|  +  a|  _ 

^A|  +  A?  +  A? 


(4.17) 


Lastly,  the  ofF-diagonal  terms  are  necessary  to  complete  the  check  for  the  Legendre-Clebsch 
condition: 


dm 

da  dp 


dm 

dp  da 


=  Ai  A  sin  a  sin  /3  —  Ay  A  cos  a  sin  /3 


(4.18) 


Using  Equations  (4.11)  and  (4.12)  to  eliminate  the  control  angles  from  the  right-hand  side 
yields  the  following: 


dm  _  dm  _  A(A^AyAi-A^AyAi) 

dadp  ~  dpda  ~  /xi  +  xyxi  +  xi  +  xj 


(4.19) 
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Thus,  the  ofF-diagonal  terms  are  identically  equal  to  zero  and  the  diagonal  terms  are  both 
greater  than  or  equal  to  zero.  This  makes  the  second  partial  derivative  matrix  positive 
semi-definite,  which  satisfies  the  Legendre-Clebsch  condition.  Thus,  we  have  the  necessary 
conditions  to  conclude  that  the  control  minimizes  the  path  integral  [7].  As  mentioned 
before,  there  is  no  case  of  interest  in  which  all  three  of  the  Lagrange  multipliers 
and  Ai  are  aU  identically  equal  to  zero  throughout  the  path  integral  over  the  transfer  time. 
Since  these  costates  premultiply  the  constraint  equations  in  the  Hamiltonian,  the  optimal 
path  would  be  independent  of  the  gravitation  and  thrust  models  if  the  velocity  costates 
were  aU  zero  throughout  the  transfer. 

The  costate  equations  are  found  using  Equation  (4.7),  and  the  results  are  as  follows: 


A.  = 

^  [Ai  (2/^  +  -  2x^)  -  3x  (Ayt/  +  Ai^)] 

(4.20) 

Ay  = 

^  [Ay  {x^  +  -  2y'^)  -  3y(XiX  +  Xiz)] 

(4.21) 

A.  = 

^  [Ai  {x^  +  y^  -  2z^)  -  3z  (XiX  +  Ayi/)] 

(4.22) 

Ai  = 

-A. 

(4.23) 

Ay  = 

—  Ay 

(4.24) 

Ai  = 

-A. 

(4.25) 

In  the  two-dimensional  problem,  the  three  final  boundary  conditions  are  the  radial 
distance  r  =  the  radial  velocity  u  =  0,  and  the  tangential  velocity  v  =  l(y/R.  This 
set  of  conditions  may  be  used  to  define  a  circular  orbit  completely  with  a  desired  direction 
of  rotation.  In  the  three-dimensional  problem,  two  more  final  conditions  are  necessary  to 
correspond  with  the  two  out-of-plane  initial  costates,  A^(0)  and  Ai(0).  If  the  three  com¬ 
ponents  of  the  angular  momentum  vector  h  are  used  along  with  the  radial  and  tangential 
velocities,  a  total  of  five  scalar  boundary  conditions  are  established.  These  quantities  will 
uniquely  determine  a  circular  orbit  with  a  desired  inclination  and  ascending  node.  Since 
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h  =  fx  the  components  of  the  angular  momentum  vector  may  be  calculated  as  follows: 


1 

II 

(4.26) 

hy  =  zx  —  xz 

(4.27) 

1 

II 

(4.28) 

The  above  equations  are  useful  for  checking  end  conditions  after  integrating  the  Cartesian 
equations  of  motion  in  three  dimensions.  If  the  final  orbit  is  circular  with  a  given  radius, 
inclination  and  ascending  node,  the  components  of  the  final  angular  momentum  may  be 
found  from  the  circular  velocity  and  angular  information. 

Now  that  the  equations  of  the  states  and  costates  have  been  derived,  they  must  be 
initialized  prior  to  numerical  integration.  The  initial  values  of  the  states  are  known,  but  the 
initial  values  of  the  Lagrange  costates  are  unknown.  However,  the  number  of  independent 
unknown  initial  costates  may  be  reduced.  Either  Aa;(0),  Ay(0),  or  A^(0)  may  be  initially 
scaled  to  unity  by  dividing  the  initial  value  of  the  Hamiltonian  by  Aa;(0),  Ay(0)  or  A5,(0). 
Although  there  is  no  mathematical  reason  to  prevent  these  initial  costates  from  being  zero 
at  the  same  time,  the  resulting  problem  would  not  be  of  practical  interest.  This  is  because 
the  resulting  optimal  trajectory  would  be  independent  of  the  initial  velocity  components, 
as  can  be  seen  from  the  Hamiltonian.  Thus,  it  is  assumed  that  A;j,(0)  ^  0,  allowing  it 
to  be  divided  through  the  Hamiltonian.  As  a  result  of  this  scaling,  Aa;(0)  =  1,  and  the 
cost  is  also  scaled  since  the  Lagrangian  becomes  1/Aa.(0),  a  new  constant.  This  has  no 
effect  on  the  final  optimal  trajectory  since  the  minimum  of  the  scaled  cost  corresponds  to 
the  minimum  time.  If  the  scaling  changes  the  sign  of  the  terms  in  the  Hamiltonian,  then 
the  optimization  problem  becomes  one  of  maximizing  the  negative  of  the  time.  This  is 
discussed  in  the  section  on  optimization  of  the  two-dimensional  problem. 

The  Lagrange  costates  associated  with  the  position  and  velocity,  A^(0)  and  Ai(0), 
are  unknown  at  the  initial  time.  However,  they  may  be  set  equal  to  zero  as  a  reasonable 
starting  guess  at  the  initial  time.  This  is  because  A^  and  A^  are  both  always  equal  to  zero 
if  the  transfer  is  confined  to  the  y  plane.  In  Table  4.1,  the  initial  values  A2(0)  and  Ai(0) 
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are  indicated  as  unknowns,  because  zero  will  not  be  the  optimal  value  in  general.  Finally, 
it  will  be  shown  in  a  later  section  that  the  initial  values  Aj^(O)  and  Ai,(0)  are  equal  if  the 
starting  orbit  is  circular.  The  three-dimensional  set  of  states  and  costates  are  initialized 
as  given  in  Table  4.1: 


Table  4.1  Initialization  for  3D  Problem 


x(0) 

=  1 

A.(0) 

=  1 

2/(0) 

=  0 

Ay(0) 

=  ? 

2(0) 

=  0 

A,(0) 

=  ? 

i(0) 

=  0 

Ai(0) 

=  A,(0) 

2/(0) 

=  1 

Ay(0) 

=  ? 

2(0) 

=  0 

Ai(0) 

=  ? 

Table  4.1  shows  a  total  of  four  unique  scalar  values  to  be  found,  and  the  optimal  flight 
time  makes  five  unknowns  for  the  three-dimensional  problem.  Five  scalar  end  conditions 
are  required  at  the  final  time  to  produce  a  square  Jacobian  matrix,  which  will  be  described 
later  in  this  chapter.  To  provide  a  total  of  five  scalar  end  conditions,  the  three  components 
of  the  final  angular  momentum  vector  are  used  in  addition  to  the  radial  and  tangential 
velocities.  If  angular  momentum  and  velocity  are  matched  correctly,  their  relationship 
guarantees  the  correct  value  of  iZ.  Thus,  R  becomes  a  redundant  end  condition,  and  does 
not  need  to  appear  explicitly.  In  fact,  any  one  of  the  six  scalar  quantities  r,  u,  t?,  hy^ 
hz  could  be  eliminated  in  this  way. 


4^2  Optimal  Control  in  Two  Dimensions 

Using  the  two-dimensional  equations  of  motion  derived  in  Chapter  3,  the  Hamiltonian 
for  this  problem  becomes: 


=  1  +  Aa,i  +  \yy  +  Xi  +  ^cosa^  +  A^  +  Asina^  (4.29) 

The  control  law  which  minimizes  the  flight  time  for  a  given  end  condition  is  found  by 
setting  dH/da  =  0,  and  leads  to: 

tana=^— (4.30) 
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Using  the  control  law  given  above,  we  obtain  the  following: 


cos  a  =  ■  = 

\/^i  +  H 

-^y 

sm  a  =  ,  = 

\/^l  +  A| 

Substituting  these  results  into  the  Legendre-Clebsch  condition  yields: 


(4.31) 

(4.32) 


d^n  ^  +  A| 

v/a|  +  A| 


>  0 


(4.33) 


Equation  (4.33)  provides  the  necessary  conditions  for  a  minimum.  It  is  interesting  to  notice 
the  similarity  between  the  above  result  and  the  equivalent  development  in  the  three  di¬ 
mensions.  In  both  cases,  the  velocity  costates  appear  in  the  denominator.  Equation  (4.33) 
could  be  further  reduced  by  division,  but  it  is  left  in  this  form  for  comparison  with  the 
three-dimensional  results. 

In  order  to  be  consistent  with  the  results  of  Bryson  and  Ho  [7],  the  Hamiltonian  may 
be  multiplied  by  —1  without  loss  of  generality.  The  Lagrangian  and  Lagrange  multipliers 
wiU  aU  change  sign.  This  leads  to  an  equivalent  control  law: 


,a„a=(^) 

We  then  have  the  following: 

cos  a  =  ■  = 

y/H  +  H 

sin  a  =  .  = 


(4.34) 


(4.35) 

(4.36) 
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When  the  Lagrangian  is  —1,  the  optimization  problem  is  one  of  maximizing  the  negative 
of  the  time  of  flight.  The  Legendre-Clebsch  condition  is  still  satisfied  as  follows: 


d^H 

da^ 


(4.37) 


which  provides  the  necessary  conditions  for  a  maximum  of  the  negative  of  the  time  of 
flight,  or  the  minimum  time  to  a  given  final  boundary  condition. 

Using  the  control  law  in  Equation  (4.34),  the  differential  equations  for  the  states  and 
costates  are  as  follows: 


(4.38) 

(4.39) 

(4.40) 

(4.41) 

(4.42) 

(4.43) 


Every  solution  of  the  differential 
time  arc  to  some  final  end  condition. 


equations  for  the  state  and  costates  is  a  minimum- 
The  states  and  costates  are  initialized  as  follows: 


Table  4.2  Initialization  for  2D  Problem 


x(0)  = 

1 

A.(0)  = 

1 

?/(0)  = 

0 

Ay(0)  = 

? 

x(0)  = 

0 

Ai(0)  = 

A,(0) 

y(o)  = 

1 

Ai(0)  = 

? 
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Table  4.2  shows  the  initialization  of  the  states  and  costates  in  the  two-dimensional 
problem.  There  are  two  unknown  initial  costates,  and  the  unknown  flight  time  makes  a 
total  of  three  scalars  to  be  found.  The  end  conditions  for  r,  ti,  and  v  are  an  appropriate 
choice  for  matching  via  the  shooting  method,  which  is  discussed  later  in  this  chapter. 

Using  the  polar  equations  of  motion  from  Chapter  3,  Equations  (3.8)  -  (3.10),  the 
Hamiltonian  is  as  follows: 

=:  1  +  XrU  +  Au  - -  +  sin  +  Xy  —  +  A  cos  (4.44) 


The  optimality  condition  yields  the  following: 

dn 


d<f> 


=  A  ( Au  cos  (t>  —  Xy  sin  ^)  =  0 


Solving  for  the  polar  thrust  angle  leads  to: 


(f)  =  tan 


-1  (K 


K\ 

k) 


Using  the  control  law  given  above,  we  obtain  the  following: 

/  Xi) 

cos  (p  = 

sin  (f)  = 


A., 


v/AlTAf 

To  find  the  costates,  we  need  the  other  result  from  Euler- Lagrange  theory: 


(4.45) 


(4.46) 


(4.47) 

(4.48) 


Ar  = 

dr  “  V  7-2  y  "  r2 

(4.49) 

A„  = 

on  _  V 

r,  —  A,-  -|-  A^ 

OU  T 

(4.50) 

II 

dH  2v  u 

X  -  Xy  1“  At; 

ov  r  r 

(4.51) 
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Table  4.3  Initialization  for  Polar  Coordinate  Problem 


r(0)  =  1 

A,(0)  =  1 

0 

II 

0 

A„(0)  =  ? 

i;(0)  =  1 

A„(0)  =  ? 

which  completes  the  set  of  optimal  control  equations  in  the  polar  case.  The  polar  states 
and  costates  are  initialized  as  shown  in  Table  4.3. 

In  the  polar  coordinate  case,  there  are  two  unknown  initial  costates  along  with  the 
unknown  final  time.  These  three  quantities  are  used  in  discussion  of  the  Jacobian  matrix 
later  in  this  chapter. 

4.2,1  Comparison  of  Cartesian  and  Polar  Hamiltonians.  At  the  initial  time 
only,  the  inertial  thrust  angle  a  may  be  used  in  the  polar  Hamiltonian  since  the  polar  and 
inertial  angles  sum  to  exactly  7r/2  at  t  =  0.  At  t  =  0,  the  spacecraft  is  located  at  the  point 
(1,0)  on  the  cc-axis.  The  spacecraft  local  horizontal  direction  is  along  the  positive  ^/-axis. 
The  inertial  angle  a  is  always  measured  to  the  thrust  vector  from  the  x-axis,  and  initially 
the  polar  angle  0(0)  is  measured  to  the  thrust  vector  from  the  y-axis.  Since  the  x  and  y 
axes  are  'k  12  radians  apart,  the  two  angles  must  sum  to  7r/2.  This  relationship  may  be 
expressed  as  follows: 


sin  0(0)  =  cos(7r/2  —  0(0))  =  cosa(O)  (4.52) 

cos  0(0)  =  sin(7r/2  —  0(0))  =  sin  a(0)  (4.53) 

The  sine  and  cosine  terms  wiU  then  be  reversed  because  of  this  relationship.  This 
substitution  allows  for  direct  comparison  of  the  terms  in  each  Hamiltonian. 

Here,  the  initial  Hamiltonian  is  expressed  in  both  sets  of  coordinates,  with  “p”  de¬ 
noting  polar  and  “c”  denoting  Cartesian: 

?fp(0)  =  1  “h  +  Ay  ^ - ^  A  cos  ^ — -  |-  A  sin  (4.54) 
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Hc{0)  =  1  +  Xxi  +  \y  +  Ax  +  ^cosa^  +  Xy  (^-^V  +  ^sina^  (4.55) 

where  all  terms  have  their  initial  values.  The  initial  conditions  for  the  polar  and  Cartesian 
cases  are: 


r(0)  = 

1 

a;(0)  = 

1 

ti(0)  = 

0  & 

2/(0)  = 

0 

t;(0)  = 

1 

i(0)  = 

0 

m  = 

1 

Using  these  initial  conditions  and  equating  the  resulting  initial  Hamiltonian  expressions 
yields: 

A^Acosa  +  XyAsina  =  A^  +  (—1  +  A  cos  a)  +  A^Asina  (4.56) 

where  all  terms  have  their  initial  values.  By  equating  coefficients  of  sina(O)  and  cosq:(0), 
the  following  three  relationships  are  obtained: 


A,(0)  = 

Ax(0) 

(4.57) 

Ai(0)  = 

A„(0) 

(4.58) 

A^(0)  = 

A.(0) 

(4.59) 

The  equality  in  Equation  (4.57),  Ay(0)  =  Ai(0),  is  the  same  relationship  shown 
in  the  initialization  Tables  4.1  and  4.2  for  the  three-  and  two-dimensional  problems.  The 
initial  conditions  used  in  each  Hamiltonian  are  expressed  in  dimensionless,  canonical  units, 
which  allow  for  great  simplification  by  substitution  of  appropriate  ones  and  zeros.  If 
physical  units  were  used  instead,  then  Ay(0)  and  Ai(0)  would  be  related  by  some  constant 
conversion  factor  based  on  the  system  of  units.  However,  this  constant  may  be  found  by 
simply  equating  the  polar  and  Cartesian  Hamiltonians  as  done  here,  only  with  the  initial 
conditions  expressed  in  the  physical  units. 

This  same  relationship  between  Ay(0)  =  Ai(0)  also  arises  in  the  minimum-fuel  prob¬ 
lem  with  impulsive  thrust,  and  may  be  derived  from  Lawden’s  primer  vector  results  for  the 
impulsive  orbital  transfer  case  [13,  18].  In  both  the  minimum-time  and  the  minimum-fuel 
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formulations,  a  starting  circular  orbit  of  unit  radius  is  required  for  the  simple  equality 
relationship  to  hold. 


4^3  Optimal  Control  under  the  Kustaanheimo-Stiefel  (KS)  Transformation 

Although  it  is  possible  to  use  the  KS  transformation  for  problems  with  three  di¬ 
mensions  [21],  only  the  two-dimensional  problem  will  be  addressed  here  for  the  purpose  of 
comparison  with  the  two-dimensional  Cartesian  case.  In  the  next  chapter,  the  initial  values 
of  the  costates  are  presented  graphically  as  functions  of  A  and  R  for  the  two-dimensional 
case.  The  three-dimensional  case  requires  two  more  parameters,  so  this  type  of  graphical 
presentation  is  not  possible.  Thus,  only  the  two-dimensional  case  under  the  KS  transfor¬ 
mation  is  necessary  for  the  graphical  comparison. 

In  the  minimum-time  problem  [7],  the  cost  functional  is  the  real  time.  As  shown 
in  Chapter  3,  the  independent  variable  is  changed  to  s  under  the  KS  transformation. 
Therefore,  the  cost  functional  becomes: 


J 


rS-Sf 

=  t=  / 

^5  =  0 


rds 


Thus,  the  Lagrangian  is  r.  The  Hamiltonian  for  this  problem  is  as  follows: 


H  —  +  ^vi 

I  2(^v)  -  p 


2{v^v)  —  p 


2r 


Ui  +  cos  7 


) 


+ 


2r 


U2 


Defining  Aj  =  At  +  1,  the  Hamiltonian  becomes: 

2{iFv)  —  p 


H  —  Auj^i  4“  Au3n2  +  A^j 


2r 


I  ] 

Uy  +  cos  7 


+  A, 


V2 


2{iFv)  —  p 


U2 


+  “  sin  7|  +  Xtr 


(4.60) 


+  sin  7I  +  r  +  Xir  (4.61) 


(4.62) 
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As  before,  the  parameter  A  is  a  function  of  time,  defined  by  A  =  T/{mo  +  rht).  T  is  the 
thrust,  mo  is  the  initial  spacecraft  mass,  and  m  is  the  constant  mass  flow  rate.  The  optimal 
control  law  is  found  by  setting  dH/d'f  equal  to  zero.  For  a  minimum,  the  result  is: 

tan  7  =  — ^  (4.63) 

This  choice  of  sign  satisfies  the  Legendre- Clebsch  condition  for  a  minimum.  From  this 
control  law,  we  have  the  following: 


cos  7 


sin  7 


(4.64) 

(4.65) 


These  relationships  may  be  used  to  eliminate  the  sin  7  and  cos  7  terms  from  the  Hamilto¬ 
nian.  The  costate  equations  are  then  found  using  the  canonical  relationship  A'  =  —dH/dq^ 
in  which  q  =  {ui^U2^Vi^V2,t).  Recall  that  the  primes  indicate  differentiation  with  respect 
to  the  fictional  time,  s.  Using  this  relationship  and  taking  the  indicated  partial  derivatives 
produces  the  following  five  first-order  differential  equations: 


a;,  = 


2(fr  u)  -  /X 


2r 


2u 


XL  = 


«3 


+  (3/2)Ar(^/^V  1/^27+^ -  2ui 

2(u^u)  —  fj, 


2r 


2U2 


+  A„,,«2)  -  A. 


V2 


+  (3/2)Ar(‘/2)u2^A2,  +  A2^  -  2u2 

A(,j  =  ( — ^T^)  +  X^^U2)  —  Au, 


(4.66) 


(4.67) 


(4.68) 

(4.69) 

(4.70) 
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If  m  =  0,  then  AJ  =  0  as  well.  In  this  case,  wiU  be  a  constant.  It  would  then  be  possible 
to  divide  the  Hamiltonian  through  by  A^  which  would  scale  the  remaining  costate  variables, 
and  eliminate  A^  from  the  problem.  If  m  /  0,  A^  must  be  retained  because  the  real  time 
wiU  explicitly  appear  in  the  equations  of  motion.  The  Hamiltonian  may  be  set  equal  to 
zero  at  the  initial  time  by  adding  an  arbitrary  constant,  since  the  constant  wiU  contribute 
nothing  to  the  partial  derivatives.  Then,  it  is  possible  to  solve  for  initial  A^  (or  At;^)  as  a 
function  of  the  initial  values  of  the  remaining  states  and  costates: 

A«(0)  =  -^^A^(0)  +  A2^(0)  (4.71) 

At  the  initial  time,  A^.^  =  2At,2-  This  may  be  shown  by  equating  the  KS  Hamiltonian 
to  another  Hamiltonian  in  polar  coordinates,  but  with  ficticious  time.  This  procedure  is 
analogous  to  the  derivation  given  earlier  using  the  Cartesian  and  polar  Hamiltonians.  Also, 
the  Hamiltonian  may  be  scaled  such  that  A^^^  =  1.  Thus,  there  are  three  remaining  values 
that  must  be  found  to  solve  the  boundary  value  problem  of  coplanar  transfer  between  two 
circular  orbits.  As  shown  in  Table  4.4,  they  are:  5/,  A„2(0),  and  Xy^{0).  If  A*  is  not  used, 
then  the  three  values  that  must  be  found  are  5/,  A„j(0),  and  Au2(0).  The  other  two  costates 
are  found  from  A,;^  =  2At,2,  and  from  solving  H{0)  =  0  for  Since  the  final  circular 
orbit  may  be  described  using  three  scalar  values,  the  number  of  unknowns  is  the  same  as 
the  number  of  end  conditions  to  be  matched. 


Table  4.4  Initialization  for  KS  Problem 


tti(O)  = 

1 

A„,(0)  =  1 

^2(0)  = 

0 

A„,(0)  =  ? 

^^i(O)  = 

0 

A,.(0)  =  2A„,(0) 

^’2(0)  = 

1/2 

A.3(0)  =  ? 

4-4  Solution  of  the  Optimal  Control  Problem 

4-4'^  The  Shooting  Method.  The  shooting  method  [16]  is  an  indirect  technique 
for  solving  a  two-point  boundary  value  problem  by  numerically  perturbing  a  reference 
trajectory.  It  is  named  after  the  classical  method  of  aiming  an  artillery  piece.  The  basic 
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idea  is  to  fire  one  sheU  as  a  reference,  then  ‘‘bracket”  the  target  by  adjusting  the  angle  up 
or  down  if  the  shot  fell  short  or  long,  respectively.  Inherent  in  this  method  is  the  need  to 
guess  the  initial  angle  for  the  first  shot,  and  the  need  for  an  effective  means  of  selecting 
the  next  angle  based  on  the  results  of  the  previous  shot. 

When  applied  to  the  continuous-thrust  spacecraft  problem,  the  differential  equations 
governing  the  state  and  costate  variables  are  numerically  integrated  to  form  the  required 
reference  trajectory.  Since  the  initial  values  of  the  costate  variables  are  usually  unknown, 
they  must  be  guessed.  This  normally  results  in  failure  to  meet  the  desired  end  conditions  for 
the  state  of  the  spacecraft,  even  though  every  trajectory  satisfies  the  optimality  condition. 
To  correct  this,  additional  trajectories  are  propagated  with  slight  changes  in  the  initial 
values  of  the  costates.  Then,  a  matrix  of  partial  derivatives  is  formed  to  quantify  the 
influence  of  initial  costate  values  on  final  states.  With  this  “secant”  or  Jacobian  matrix, 
the  optimal  initial  values  of  the  costates  may  be  found  approximately.  If  the  elements  of 
the  Jacobian  matrix  are  found  numerically  with  one-sided  difference  approximations,  as 
is  done  in  this  implementation,  the  search  is  a  gt/asz-Newton  method.  Also,  the  Jacobian 
matrix  represents  only  a  first  order  linear  approximation  for  a  typically  highly  nonlinear 
set  of  equations.  Often,  the  nonlinearity  results  in  a  process  that  will  not  converge  from 
poor  initial  guesses.  The  Jacobian  matrix  may  also  be  found  analytically,  allowing  for 
integration  of  the  equations  of  variation  and  better  convergence  properties.  This  technique 
is  a  multi- variable  application  of  Newton’s  method  [14,  16]  to  a  nonlinear  problem,  which  is 
only  guaranteed  to  converge  within  a  “small”  neighborhood  of  the  solution  point.  The  size 
of  this  neighborhood  is  determined  by  the  extent  over  which  the  problem  is  approximately 
linear  about  the  solution  point.  Since  the  exact  equations  for  the  states  and  costates  do 
not  have  an  analytical  solution,  the  only  way  to  approximate  the  size  of  the  nearly  linear 
region  is  through  numerical  methods.  Appendix  A  contains  a  flowchart  that  outlines  the 
shooting  method  algorithm.  Appendix  A  also  contains  a  description  of  the  convergence 
criteria  used  for  aU  of  the  numerical  examples  in  this  dissertation. 

When  attempting  to  solve  a  problem  of  this  nature,  the  analyst  has  a  choice  of 
variables  to  use  in  forming  the  Jacobian  matrix.  It  is  generally  a  good  idea  to  keep  the 
partial  derivative  matrix  square,  since  the  inverse  will  be  required  to  find  the  corrections 


4-18 


to  the  initial  costates.  If  the  Jacobian  matrix  is  not  square,  other  alternatives  exist  such 
as  left  and  right  inverses  [22].  If  the  final  time  of  a  boundary  value  problem  is  unknown, 
but  the  final  states  are  specified,  then  the  final  time  may  be  chosen  as  an  input  to  the 
Jacobian  matrix.  This  way,  the  unknown  final  time  may  be  refined  from  an  initial  guess 
along  with  the  costates. 

The  first  step  in  solving  this  boundary  value  problem  is  to  guess  the  two  unknown 
costates  and  a  final  time,  then  numerically  integrate.  Once  the  final  time  has  been  reached, 
the  final  state  values  are  examined  to  see  how  close  they  came  to  the  desired  final  conditions. 
Since  they  will  undoubtedly  be  wrong,  the  initial  costates  and  final  time  must  be  adjusted 
to  try  again.  If  this  process  is  convergent,  the  final  conditions  wiU  be  met  to  a  desired 
accuracy  after  several  iterations.  Time  is  the  independent  variable,  so  a  final  value  of  the 
time  makes  a  reliable  stopping  condition  for  integrating  the  equations  of  motion.  The 
stopping  condition  r  =  i?  is  not  as  useful,  since  the  equations  of  motion  apply  equally  to 
orbit-lowering  and  orbit-raising.  Thus,  some  initial  values  of  the  costates  lead  to  a  decrease 
in  r,  so  that  r  =  R  may  not  occur  at  or  before  t  =  tj.  Therefore,  refining  the  final  time  is 
the  most  practical  approach  to  the  problem. 

To  adjust  the  initial  costate  values  and  final  time,  we  form  the  Jacobian  matrix  of 
partial  derivatives  numerically  to  see  how  the  final  state  errors  depend  on  initial  guesses. 
We  will  now  describe  the  formulation  of  the  quasi-Newton  step  using  the  polar,  two- 
dimensional  case  which  results  in  a  3  X  3  Jacobian  matrix  since  there  are  two  initial  costates 
and  the  time  of  flight  to  be  determined.  By  comparison,  the  three-dimensional  case  wiU 
have  a  5  X  5  Jacobian  matrix  corresponding  to  four  unknown  initial  costates  and  the  time 
of  flight. 

If  the  initial  value  of  A,,  was  not  taken  to  be  unity  as  explained  before,  a  reasonable 
formulation  for  the  quasi-Newton  step  in  polar  coordinates  would  be  as  follows: 

AA,.(0)  dr{tj)/dXr{0)  dT{tf)/dX^{0)  dr{tf)/dXy{0)  Ar(t/) 

AA,,(0)  =  du{tf)/dXr{0)  du{tf)/dX^{0)  du{tf)/dX^{0)  Au{tf)  (4.72) 

AA^(O)  dv{tf)/dXr{0)  dv{tj)/dX^{0)  dv{tj)/dXy{0)  Au(t/) 
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Each  of  the  partial  derivatives  is  formed  by  making  small  changes  to  the  initial  conditions 
of  a  reference  trajectory,  and  noting  the  changes  in  the  final  conditions.  As  mentioned 
previously,  the  partial  derivatives  are  approximated  numerically,  so  this  technique  is  a 
qwosi-Newton  method. 

Since  the  initial  value  A,.(0)  may  be  scaled  to  unity,  there  is  no  need  to  include  it  in 
the  above  formulation.  However,  the  final  time  is  also  unknown  in  this  problem.  Rather 
than  make  manual  changes  in  final  time,  it  is  possible  to  incorporate  it  directly  into  the 
quasi'Newton  step  formulation: 


Atf 

dT(tj)/dtj  dTit;)/dX^{0)  dT{t,)/dX,{0)' 

-X 

Ar{tf) 

AA„(0) 

= 

duitj)/dt^  du{t;)/dX^{0)  du{t,)/dX^iO) 

Au{tf) 

AA„(0) 

dv{tf)/dtj  dv{tj)/dx^(0)  dv{tf)/dx,i0) 

Av{tj)  _ 

This  is  an  original  formulation  that  saves  computing  time,  since  it  completely  auto¬ 
mates  the  shooting  technique.  In  this  form,  the  boundary  value  problem  may  be  imbedded 
in  a  larger  programming  loop  that  allows  for  parameter  variations.  This  way,  we  can 
more  easily  study  the  effects  of  changing  the  thrust  level,  mass  flow  rate,  final  radius,  and 
gravitational  constant. 

The  above  Jacobian  matrix  formulation  automates  the  search  over  a  range  of  fixed 
final  time  values.  Each  “shot”  in  the  shooting  method  is  still  a  minimum-time  arc,  but 
this  technique  allows  us  to  match  the  desired  end  conditions  and  determine  the  correct 
final  time  simultaneously.  In  the  three-dimensional  case,  the  quasi- Newton  step  becomes: 


Atf 

Au{tf) 

AAj,(0) 

Av(tf) 

AA,(0) 

Ahg(tf) 

AAy(O) 

Ak,(tf) 

AAi(O) 

(4.74) 
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The  Jacobian  matrix  is  denoted  by  J.  This  expression  corresponds  to  the  results  given 
in  Table  4.1.  The  partial  derivatives  in  the  5x5  Jacobian  matrix  are  similar  to  the 
two-dimensional  case,  using  the  variables  given  in  the  above  equation. 


Dynamic  Step  Limiter,  To  improve  the  convergence  properties  for  these 
problems,  an  original  modification  is  made  to  Newton’s  method  by  providing  a  variable 
scaling  factor  for  the  quasi-Newton  step.  Let  P  be  a  vector  of  unknowns  such  as  P  = 
(t/,  Au(0),  A^(0)),  and  let  AP  be  the  quasi-Newton  step  from  Equation  (4.73),  where 


AP  = 


At  j 

AA,(0) 

AA,(0) 


Then,  the  following  equation  is  used: 


^t+i  —  Pi  + 


AP, 


1  + 


AP 


(4.75) 


(4.76) 


This  modification  provides  a  dynamic  scaling  effect  on  the  quasi-Newton  step.  If 
is  large,  the  scaling  limits  the  individual  component  changes  to  less  than  unity,  without 
changing  the  step  direction.  If  Ap  is  small,  the  scaling  does  not  affect  the  magnitudes  of 
the  individual  component  changes  very  much,  so  the  quadratic  convergence  rate  is  nearly 
maintained  in  the  neighborhood  of  the  solution  point.  The  scaling  preserves  the  direction 
of  the  quasi-Newton  step  in  the  search  space,  but  prevents  the  magnitude  from  becoming 
unreasonably  large,  as  can  happen  with  nonlinear  problems. 


A  proof  of  convergence  for  Newton’s  method  is  given  by  Luenberger  [14]  for  the  mul¬ 
tivariable,  nonlinear  case.  To  quote  the  source,  ‘Hhe  theorem  can  be  paraphrased  roughly 
by  simply  saying  that  Newton’s  method  converges  provided  that  the  initial  approximation 
Xi  is  sufliiciently  close  to  the  solution  Xq.”  Luenberger  also  states,  “One  device  useful  in 
these  situations  is  to  begin  iterating  with  a  slower  but  surer  technique  and  then  change 
over  to  Newton’s  method  to  gain  the  advantage  of  quadratic  convergence  near  the  end  of 
the  process.”  This  is  exactly  the  purpose  of  the  dynamic  scaling  technique  presented  here. 
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From  numerical  experience,  it  has  been  found  that  an  initial  approximation  that  is  just 
outside  of  the  region  of  convergence  will  often  result  in  a  very  large  quasi-Newton  step 
magnitude.  By  limiting  the  magnitude  of  the  step,  it  is  possible  to  exploit  the  directional 
information  from  the  quasi-Newton  step  while  potentially  increasing  the  region  of  conver¬ 
gence.  Due  to  the  nonbnearity  of  the  problem,  this  result  has  only  been  shown  to  hold  true 
to  varying  degrees  based  on  nunierical  experimentation.  The  one  definite  statement  that 
can  be  made  about  this  scaling  technique  is  that  it  cannot  make  the  region  of  convergence 
smaller.  If  the  quasi-Newton  step  is  denoted  by  AP  as  before,  the  convergence  proof  [14] 
requires  that: 

<77  (4.77) 

where  77  is  some  sufficiently  small  scalar  quantity  to  provide  convergence.  Clearly,  it  may 
be  seen  that 

<  APi  <r)  (4.78) 

The  modified  Newton  step  satisfies  the  proof,  as  long  as  the  original  Newton  step  does. 
The  modified  quasi-Newton  method  may  take  more  iterations  to  acheive  convergence,  but 
it  also  increases  the  radius  of  convergence,  based  on  numerical  experience. 

In  this  research,  the  modified  quasi-Newton  method  is  used  for  its  simplicity.  The 
initial  costate  approximations  contain  an  entire  control  history  in  a  compact  form,  since  the 
thrust  angle  is  produced  automatically  with  numerical  integration  through  the  final  time. 
Thus,  by  applying  the  modified  quasi-Newton  method  to  determine  the  initial  costates 
and  time  of  flight,  it  is  possible  to  explore  a  large  range  of  control  functions.  Further 
extensions  of  Newton’s  method  are  available  [9],  but  they  add  more  complication  to  the 
search  technique.  Most  importantly,  the  purpose  of  this  research  is  to  provide  approximate 
models  with  enough  accuracy  so  that  the  modified  quasi-Newton  method  converges  without 
further  modification.  A  convergence  sensitivity  study  is  presented  in  Chapter  5,  which 
shows  the  number  of  iterations  required  for  convergence  for  the  entire  practical  range  of  R 
and  A. 

Figure  4.1  shows  a  comparison  of  the  modified  quasi-Newton  method  (MNM)  with 
the  unmodified  quasi-Newton  method  (NM).  This  is  a  a  well-known  Earth-to-Mars  transfer 
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Modified  quasi-Newton  vs.  quasi-Newton  Method 
for  the  Bryson-Ho  Example 


Figure  4.1  Iteration  History  for  Modified  vs.  Unmodified  Quasi-Newton  Method 


example  developed  by  Bryson  and  Ho  [7],  which  will  be  discussed  in  greater  detail  in 
Chapter  6.  The  iteration  history  shows  seven  modified  quasi-Newton  steps  with  solid 
lines,  and  six  unmodified  quasi-Newton  steps  with  dashed  lines.  The  solid  line  segments  in 
Figure  4.1  are  seen  to  have  slope  magnitudes  less  than  or  equal  to  the  dashed  line  segments 
between  iteration  numbers  1  and  2.  Since  iterations  2  through  6  start  from  different  values 
using  the  two  methods,  the  modified  quasi-Newton  steps  may  have  slopes  that  are  larger, 
smaller  or  the  same  as  the  quasi-Newton  steps.  In  this  example,  the  modification  caused 
an  increase  of  one  iteration  over  the  quasi-Newton  method  to  reach  the  same  converged 
values.  Appendix  A  describes  the  numerical  criteria  used  in  this  and  all  other  examples 
presented  in  this  dissertation. 

Table  4.5  shows  a  case  with  R  =  2.2  and  A  =  0.01  in  which  the  quasi-Newton 
method  begins  to  diverge  on  the  third  iteration.  The  modified  quasi-Newton  method 
achieves  convergence  in  11  iterations. 


4^5  Summary 

This  chapter  presents  the  minimum-time  optimal  control  formulation  in  several  coor¬ 
dinate  systems  for  the  orbital  transfer  problem  considering  gravity  and  continuous  thrust. 
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Table  4.5  Comparison  of  Search  Methods  for  R  =  2.2,  A  =  0.01 


quasi-Newton  Method  Modified  quasi-Newton  Method 


it. 

A„(0) 

A.(0) 

time 

it. 

A„(0) 

A„(0) 

time 

1 

0 

1 

32.58001 

1 

0 

1 

32.58001 

2 

-0.01441 

0.637370 

32.44273 

2 

-0.00295 

0.925693 

32.55188 

3 

2.19319 

0.912565 

59.36023 

3 

-0.00138 

0.825579 

32.92815 

4 

-51.25330 

95.665720 

-54.05710 

4 

0.00441 

0.823785 

33.02219 

divergence 

5 

0.033265 

0.816452 

33.48742 

6 

0.059477 

0.817559 

33.90471 

7 

0.080149 

0.823066 

34.24611 

8 

0.093224 

0.829128 

34.47590 

9 

0.098446 

0.832526 

34.57503 

10 

0.099224 

0.833181 

34.59116 

11 

0.099240 

0.833198 

34.59154 

convergence  achieved 

In  an  original  presentation,  the  minimum-time  formulation  is  also  applied  to  the  equations 
of  motion  under  the  KS  transformation.  Numerical  initialization  techniques  are  discussed 
for  each  set  of  optimal  control  equations.  The  shooting  method  is  described  as  a  means 
to  solve  the  boundary  value  problem.  Finally,  an  original  dynamic  scaling  modification  to 
the  quasi-Newton  method  is  provided  which  improves  the  convergence  properties  of  the 
shooting  method.  An  example  is  provided  showing  differences  in  the  convergence  proper¬ 
ties  of  the  quasi-Newton  method  and  the  modified  quasi-Newton  method.  The  modified 
quasi-Newton  step  satisfies  the  conditions  for  convergence,  as  long  as  the  conditions  are 
met  by  the  unmodified  quasi-Newton  step,  and  often  when  the  unmodified  quasi-Newton 
step  does  not. 
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V.  Optimal  Initial  Costate  Locus 


In  this  chapter,  the  optimal  initial  Lagrange  multipliers  are  modeled  as  functions 
of  the  problem  parameters  which  include  R  =  the  final  radius,  and  A  =  A(0),  the 

initial  thrust  acceleration.  This  is  accomplished  by  first  examining  the  functional  form 
of  the  costates  graphically,  then  dividing  the  resulting  locus  into  three  distinct  regions 
for  separate  analysis.  A  combination  of  analytical  and  empirical  techniques  is  used  to 
model  the  regions  of  the  costate  locus.  The  models  are  then  evaluated  by  measuring 
the  convergence  sensitivity  for  the  entire  practical  range  of  A  and  R.  Finally,  the  initial 
costate  locus  is  presented  under  the  KS  transformation  to  provide  a  qualitative,  graphical 
comparison  to  the  Cartesian  form.  Appendix  A  contains  a  description  of  the  numerical 
convergence  criteria  used  to  determine  the  “exact”  solutions  that  appear  in  the  figures  and 
examples. 

To  illustrate  the  importance  of  developing  the  initial  costate  models,  one  may  simply 
consider  the  alternatives.  There  are  no  models  available  in  the  literature  that  provide 
initial  costate  estimates  for  the  minimum-time,  continuous  thrust  orbit  transfer  problem 
as  functions  of  the  problem  parameters.  As  mentioned  previously,  others  [15, 18]  have  tried, 
with  limited  success,  to  use  the  Hohmann  transfer  to  initialize  the  Lagrange  costates  for 
the  minimum-fuel  problem  with  coasting  arcs.  However,  the  minimum-fuel  problem  allows 
for  throttling  and  has  a  Lagrangian  based  on  fuel  mass,  neither  of  which  appear  in  the 
minimum-time  problem.  Thus,  the  minimum-fuel  results  provide  no  reliable  information 
to  initialize  the  costates  in  the  minimum-time  problem. 

The  only  way  to  solve  the  minimum- time  boundary  value  problem  is  to  guess  some 
initial  Lagrange  costate  values  and  hope  for  good  fortune,  resulting  in  a  convergent  case. 
However,  good  fortune  is  notoriously  unreliable.  Even  if  one  happens  upon  a  convergent 
case  to  the  desired  final  radius,  it  is  unlikely  that  the  parameter  value  of  the  thrust  accel¬ 
eration  will  match  the  spacecraft  design.  Then,  the  task  is  to  adjust  the  value  of  A  until  it 
matches  the  desired  value.  If  A  is  too  large,  it  may  be  reduced  by  some  small  percentage, 
and  used  with  the  initial  costate  values  from  the  known  case.  If  the  reduction  percentage 
is  small  enough,  the  problem  may  converge  for  the  new  value  of  A.  This  process  may  be 
repeated  until  the  desired  value  of  A  is  achieved.  The  drawback  to  this  technique,  known 
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as  the  continuation  method  [12],  is  that  the  step  size  for  A  to  result  in  convergence  is  not 
constant,  and  one  must  therefore  “hand-hold”  the  process  by  monitoring  the  convergence 
behavior  through  perhaps  hundreds  of  small  changes  in  A.  Depending  on  the  available 
computing  power  and  sensitivity  of  the  problem,  this  process  could  take  days  or  weeks  to 
complete.  Again,  one  must  also  have  found  a  “seed”  case  to  begin  the  search. 

Suppose  =  2,  and  we  have  a  converged  case  for  A  =  100.  A  numerical  description 
of  convergence  for  this  and  aU  other  examples  is  presented  in  Appendix  A.  If  A  is  multiplied 
by  0.9,  then  convergence  may  be  achieved  using  the  last  values  of  the  initial  costates  as 
the  new  starting  guess.  If  this  procedure  is  repeated  more  than  roughly  20  times,  the 
shooting  method  will  take  more  iterations  to  converge,  and  eventually  will  not  converge  at 
aU.  Then,  the  multiplication  factor  must  be  increased  to  perhaps  0.95,  and  higher  still  as 
A  decreases.  This  phenomenon  is  due  to  the  sensitivity  of  the  system  to  initial  conditions, 
which  increases  for  the  increased  flight  times  associated  with  small  thrust  values  at  a 
given  i2.  Although  this  process  is  difficult  and  time  consuming,  it  does  produce  valuable 
information.  After  completing  the  process  for  a  large  range  of  A,  it  is  instructive  to 
make  a  plot  of  the  converged  initial  values  of  the  costates  versus  one  another  to  examine 
their  behavior.  This  is  an  original  presentation  technique  [23,  24]  which  first  appears  in 
Figure  5.1.  This  choice  of  axes  is  motivated  by  the  polar  thrust  angle  (f)  as  shown  in 
Figure  3.2.  The  initial  value  of  ([)  may  be  measured  from  the  Av(0)  axis  to  a  point  on  the 
locus  with  the  vertex  at  the  origin.  Thus,  the  initial  thrust  angle  may  be  seen  relative  to 
the  spacecraft  local  horizon  directly  from  the  figure.  For  large  values  of  A,  the  initial  thrust 
angle  ^(0)  approaches  90  degrees,  and  for  small  values  of  A,  ^(0)  approaches  zero  degrees, 
in  the  spacecraft  horizontal  direction.  The  analysis  of  the  initial  costate  locus  is  carried 
out  using  the  two-dimensional  Cartesian  coordinate  formulation,  so  it  is  important  to 
remember  that  Ai(0)  =  Au(0)  and  Ay(0)  =  Ai,(0),  as  shown  in  Equations  (4.58)  and  (4.59). 

As  the  thrust  acceleration  decreases,  the  optimal  trajectory  will  cover  a  larger  transfer 
angle  to  reach  the  same  final  radius.  There  are  particular  values  of  A  that  result  in 
integer  values  of  orbital  revolutions,  and  the  points  corresponding  to  one  through  twelve 
revolutions  are  indicated  in  Figure  5.2  for  =  2.  The  first  point  marked  with  a  circle  on  the 
outermost  curve  of  the  spiral  corresponds  to  a  one-revolution  transfer.  The  next  indicated 
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Figure  5.3  Optimal  Initial  Costate  Locus  for  =  2,  m  =  —0.01,— 0.1 

point  inward  towards  the  center  of  the  spiral  corresponds  to  a  two-revolution  transfer,  and 
so  on,  until  the  last  indicated  point  closest  to  the  center  of  the  spiral  represents  a  twelve- 
revolution  transfer.  The  number  of  revolutions  will  continue  to  increase  as  the  thrust  is 
reduced,  but  the  flight  time  increases  as  weU.  This  pattern  continues  indefinitely  if  the 
mass  flow  rate  is  equal  to  zero,  because  the  spacecraft  will  never  run  out  of  propellant 
mass.  The  values  of  A  that  correspond  to  integer  revolution  transfers  wiU  change  with  R, 
but  they  have  been  found  to  give  a  negative  value  of  <^(0)  in  aU  cases. 

The  process  used  to  create  the  solid  locus  in  Figure  5.1  is  repeated  for  different 
values  of  the  mass  flow  rate,  and  shown  in  Figure  5.3.  When  rh  =  —0.1  mass  units  per 
time  unit,  the  final  time  cannot  exceed  10  time  units,  since  aU  of  the  spacecraft  mass 
would  be  consumed.  Similarly,  the  final  time  cannot  exceed  100  time  units  for  m  =  —0.01. 
Thus,  the  locus  spiral  will  stop  for  some  minimum  value  of  A,  because  small  values  of  A 
correspond  to  long  flight  times.  When  rh  is  taken  to  be  zero,  there  is  no  time  limit  for 
the  powered  trajectory,  and  any  value  of  A  may  be  used.  Therefore,  the  m  =  0  locus 
may  continue  spiraling  indefinitely  as  /I  is  reduced.  It  is  clear  from  Figure  5.3  that  the 
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changes  in  m  have  little  effect  on  the  location  of  the  initial  costates.  A  possible  physical 
interpretation  of  this  relative  insensitivity  to  m  is  that  mass  flow  rate  has  no  immediate 
effect  at  the  initial  time,  since  the  spacecraft  mass  is  normalized  to  unity  at  ^  =  0.  Also, 
only  small  changes  in  the  costates  at  the  initial  time  are  required  for  significant  changes  in 
the  trajectory  at  the  final  time,  due  to  the  sensitivity  of  the  problem.  Based  on  observations 
from  the  numerically  generated  initial  costate  loci,  the  quantity  m  is  assumed  to  be  zero  in 
the  development  of  the  initial  costate  models,  in  order  to  simplify  the  analysis.  However, 
the  equations  of  motion  always  include  m,  so  the  examples  in  Chapter  6  may  be  given 
realistic  values  of  the  mass  flow  rate. 

The  next  step  in  the  analysis  is  to  take  the  numerical  results  of  many  converged  cases 
with  R  and  A  as  parameters,  and  plot  the  optimal  initial  costates  as  functions  of  R  and  A, 
forming  the  loci  shown  in  Figure  5.4.  In  this  way,  the  functional  behavior  of  the  costates 
is  easily  seen  to  be  represented  by  three  distinct  regions.  Near  the  origin  of  the  A,^(0), 
A«(0)  plane,  the  optimal  initial  costates  lie  on  a  nearly  parabolic  arc.  As  A  decreases  or  R 
increases,  the  locus  moves  away  from  the  origin  on  a  nearly  elliptical  path,  and  eventually 
spirals  into  the  point  =  1,  =  0.  This  point  represents  the  limiting  case  of  purely 

tangential  thrust  [26].  The  optimal  initial  costate  loci  show  a  common  tendency  to  spiral 
towards  the  point  (1,0)  with  m  =  0.  Two  example  points  are  shown  in  Figure  5.4,  which 
correspond  to  an  anti-satellite  avoidance  mission  (AS AT)  [10]  and  the  well-known  Earth- 
to-Mars  transfer  example  given  by  Bryson  and  Ho  [7].  These  examples  will  be  described 
in  detail  in  Chapter  6.  The  parameter  S  =  y/{R  —  1)/ A  is  developed  in  the  next  section. 

The  parabolic,  elliptic,  and  spiral  regions  of  the  costate  locus  will  be  addressed 
separately  for  modeling  purposes.  The  parabolic  region  represents  orbital  transfers  that 
take  less  than  about  one  quarter  revolution  to  complete,  either  due  to  high  thrust  or  small 
radius  change.  This  situation  lends  itself  to  an  analytical  approach,  which  is  presented  in 
the  next  section. 

5.1  Parabolic  Region 

For  orbital  transfers  with  large  A  or  small  72,  the  optimal  trajectory  tends  to  be  a 
nearly  linear  path  which  takes  less  than  roughly  one  quarter  revolution,  based  on  numerical 
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Figure  5.4  Optimal  Initial  Costate  Loci 

investigation.  Under  these  conditions,  the  problem  may  be  approximated  with  gravity-free 
space,  since  it  is  the  influence  of  gravity  that  causes  multiple  revolution  trajectories. 

Large  A  and  small  R  values  correspond  to  the  parabolic  region  of  Figure  5.4,  near 
the  origin.  If  the  equations  of  motion  are  approximated  by  setting  //  =  0  and  m  =  0,  the 
differential  equations  for  the  states  and  costates  may  be  integrated  in  closed  form.  Although 
ignoring  the  gravity  may  not  seem  to  be  the  most  obvious  choice  for  approximation,  it  does 
reduce  the  problem  to  a  system  of  algebraic  equations.  Further,  the  boundary  conditions 
are  stiU  chosen  to  be  circular  orbits.  Thus,  gravity  stiU  has  an  influence  on  the  solution 
since  the  boundary  conditions  depend  on  the  nominal  gravitational  constant  value,  fi  =  1. 

Once  these  two  approximations  (//  =  0,  m  =  0)  have  been  made,  the  equations  of 
motion  in  Cartesian  coordinates  simplify  enough  to  allow  analytical  integration.  The  result 
is  a  system  of  eight  algebraic  equations  for  the  position  and  velocity  components  and  their 
associated  costates.  These  solutions  are  functions  of  A,  R,  time  and  eight  constants  of 
integration.  Recall  that  the  polar  and  Cartesian  optimal  initial  costates  are  related  by 
Equations  (4.58)  and  (4.59). 
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5AA  Equations  of  Motion  with  Zero  Gravity.  Setting  the  gravitational  constant 
to  zero,  the  differential  equations  of  motion  simplify  to: 


y  = 

— 

Ay  = 

A-p  — 

Ay  = 

The  last  four  equations  may  be  integrated  immediately  and  substituted  into  the  first  four 
equations.  Defining  Xpo,  =  and  A„e/  =  \Jxi  +  A|,  the  fi  =  0  solutions  for  the 

costate  and  state  equations  are: 


A,  =  a  (5.7) 

A,  =  b  (5.8) 

Xx  =  —at  +  c  (5-9) 

A^  =  -bt  +  d  (5.10) 


a;  =  ^  c  -  26^c  +  3abd  -  aA^„  t) 

+  3a6^c^  +  Aa^bcd  —  2b^cd  —  a^d^  +  2aPd^  +  2Xp^,b  {be  —  ad)  ij 

X  In  ^—ac  —  bd  A|^^t  +  A^o^A  vel  )]  •\-kit  +  fca  (5.11) 


if 

1  (5-1) 

+  ■'Ij 

.  _) 

(5.2) 

0 

(5.3) 

0 

(5.4) 

-A. 

(5.5) 

-Av 

(5.6) 
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y 


=  A 


X 


X  — 


y  = 


^^pos 

+  — ^ —  i^—Za?bd?  +  Aab^cd  —  2a^cd  —  b^c^  +  2a^bc^  +  (ad  —  6c) 

In  ac  —  bd  A  ^pos^  "t"  +  ^4 


(5.12) 


A|„, 


^^pos  ^vel  +  b{bc-ad)lii(^-ac-bd+Xl^,t  +  Xpo,Xvei)]  +  ki  (5.13) 
-bXpa,X^ei  +  a  {ad  -  6c)  In  [-ac  -bdA  X^^^t  +  Ap<„A„e/)]  +  ^2  (5.14) 


The  terms  a,  6,  c,  d,  Ajj,  k2^  and  k^  are  constants  of  integration.  Because  the  Lagrange 
multipliers  appear  in  the  Hamiltonian  as  linear  terms,  the  initial  value  of  one  of  them  may 
be  scaled  to  unity.  For  this  system,  we  choose  Aa;(0)  =  1,  thus  a  =  1.  Also,  if  the  initial 
state  is  on  a  circular  orbit,  we  have  b  —  diS  shown  in  Chapter  4  by  equating  the  system 
Hamiltonian  expressed  in  polar  and  Cartesian  coordinates. 

The  constants  ki  and  k2  may  be  eliminated  from  Equations  (5.13)  and  (5.14)  by 
using  the  velocity  component  end  conditions,  i(0),  i(^/),  ^(0),  y(tj)  and  the  final  time,  tj. 
The  final  velocity  components  come  from  the  desired  final  orbit.  It  should  be  noted  that 
Xpos  does  not  have  a  time  argument  since  it  is  a  constant  in  the  zero-gravity  case  where 
Xpos  =  a/1  +  6^,  using  the  definition  given  previously.  Performing  these  operations  yields: 


(x(<j)  x(0))A°.,  ^  A,„(A.„(0)-A..,(t,))  +  i(d-6’)i  (6.16) 

=  V,(A..,(0)-A..,(,;))+(^)i  (6.16) 


L 


(  —b  —  bd+  Apo,A„e)(0)  \ 

y  -6  —  6d  +  Ap,,,  Atie/(0)  +  Xj^^tf  j 


(5.17) 


5.1.2  Rectilinear  Case.  In  the  special  case  in  which  the  initial  and  final  velocity 
components  are  all  equal  to  zero,  the  optimal  trajectory  is  a  straight  line.  This  rectilinear 
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case  allows  for  a  simplification  of  the  above  equations.  They  take  the  form: 

0  =  A,,.  (A„,(0)  -  XUW)  +  b{d-b^)L  (5.18) 

0  =  A,„  (A.,,(0)  -  A„*,(f^))  +  L  (6.19) 

The  only  difference  in  the  above  two  equations  is  the  coefficient  of  the  logarithmic  term, 
so  these  coefficients  must  be  equal.  The  logarithmic  term  is  not  zero  unless  the  final  time 
is  zero,  but  this  is  a  degenerate  case  since  there  would  be  no  transfer  at  all.  Also,  the  term 
Apoj  is  greater  than  or  equal  to  one,  since  Aj,<„  =  v^TTP^.  The  equality  between  coefficients 
leads  to  the  following  relationship: 

(6*  +  l)  (d-fc^)  =0  (5.20) 

Clearly,  the  real  solution  is  d  =  b^.  If  this  result  is  substituted  back  into  Equation  (5.18) 
or  (5.19),  the  logarithmic  term  vanishes.  Since  A,,,  >  1  by  definition,  we  have: 


Av*((0)  =  A.,,(f;)  (5.21) 

Using  the  definition  A„,i  =  y'^A|  +  A^  and  the  solutions  for  A^  and  Ay  given  in  Equar 
tions  (5.9)  and  (5.10),  this  becomes: 

Vb^  +  cP  =  +  {d-btjf  (5.22) 

Squaring  both  sides,  expanding,  collecting  terms  and  using  d  =  ft*  leads  to: 

0  =  _2ft^ +t;ft*  -  2ft  +  r,  (5.23) 

0  =  (ft*  +  1)  (ft  -  i^/2)  (5.24) 

Again  taking  the  real  solution,  the  result  is  ft  =  t//2.  As  mentioned  previously,  Euler- 
Lagrange  theory  does  not  provide  enough  information  to  solve  the  orbital  transfer  boundary 
value  problem.  This  situation  does  not  change  by  approximating  gravity  to  be  zero.  Thus, 
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additional  information  is  required  to  solve  the  algebraic  equations.  Newton’s  law  may  be 
used  to  approximate  the  optimal  time  of  flight,  allowing  for  a  solution.  Having  assumed 
m  =  0,  the  thrust  will  switch  directions  midway  through  the  trajectory  to  decelerate  to  a 
stop,  and  the  time  of  flight  is  tj  =  We  now  define  an  original  parameter, 

5,  as  follows: 

5  =  ^{R-l)/A  (5.25) 

The  parameter  S  has  the  canonical  units  of  time,  and  is  equal  to  one-half  of  the  flight  time 
on  a  straight-line  (rectilinear)  trajectory  in  field  free  space  with  stationary  end  conditions. 
Since  it  includes  both  R  and  A,  the  S  parameter  provides  a  convenient  way  to  quantify 
the  regions  of  the  costate  locus,  as  shown  in  Figure  5.4.  Since  Ai(0)  =  b  and  A^(0)  =  d, 
we  have  an  original  solution  for  the  rectilinear  case: 

Xi{0)  =  S  (5.26) 

A^(0)  =  (5.27) 

if  =  2S  (5.28) 

These  solutions  are  for  the  simplest  case  of  no  gravity,  no  mass-flow  rate,  and  zero  velocity 

end  conditions.  Note  that  A^  =  bo  and  Ay  =  do  define  a  parabola  in  the  Ay,  A^,  plane,  or 

equivalently  in  the  A^,  A^  plane  as  shown  by  Equations  (4.58)  and  (4.59). 

Equations  (5.26),  (5.27),  and  (5.28)  along  with  b  =  c  and  a  =  1  provide  an  approxi¬ 
mate  analytical  solution  for  the  initial  values  of  the  Lagrange  costates  and  the  final  time. 
These  approximations  can  be  used  to  start  sub-optimal  transfers,  or  as  starting  points  for 
solving  the  two-point  boundary  value  problem.  Further  refinement  is  possible  by  scaling 
these  results,  as  will  be  shown  next. 

In  Figure  5.5,  the  exact  initial  costates,  obtained  numerically  as  described  in  Ap¬ 
pendix  A,  are  shown  on  the  solid  line  for  the  zero-gravity  case.  These  are  obtained  by 
numerically  solving  the  fi  =  0  =  m  two  point  boundary  value  problem  with  circular  end 
conditions  corresponding  to  an  Earth-Mars  transfer  (R=L525).  The  o  symbols  correspond 
to  A  values  of  1000,  100, 10  and  1.6,  with  the  largest  values  near  the  origin.  The  locus  data 
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Figure  5.5  Comparison  Costate  Solutions  to  Exact  Cases 

is  stopped  at  ^  =  1.6,  because  it  has  been  numerically  determined  that  the  R  =  1.525, 
fi  =  0  locus  no  longer  resembles  the  fx  =  1  locus  for  A  <  1.6.  The  dashed  line  is  a  parabola 
defined  by  Equations  (5.26)  and  (5.27)  for  bo  =  Ai(0)  and  do  =  Ay(0).  The  X  symbols 
correspond  to  the  same  values  of  A  as  the  o  symbols.  Note  the  excellent  agreement  for 
large  A,  as  expected. 

In  Figure  5.5,  the  initial  costate  approximations  based  on  Equations  (5.26)  and  (5.27) 
lie  on  nearly  the  same  parabolic  arc  as  the  exact  values  for  //  =  0,  except  that  they  appear 
to  be  shifted  upwards  along  the  arc.  This  makes  them  too  large  for  the  values  of  A,  so 
if  they  are  scaled  down  along  the  parabolic  arc,  they  wiU  more  closely  match  the  correct 
values.  An  original  scaling  approach  leads  to  improved  agreement  for  the  lower  values  of 
A.  The  scaled  points  are  shown  with  the  +  symbols,  and  the  scaling  factor  q  is  defined  as 
follows: 


5A1 


b  ^  qS 


(5.29) 


d  «  (5.30) 

It  is  interesting  to  observe  that  the  scaled  approximation  for  A  =  1.6  is  very  close  to  the 
exact  case  with  //  =  1.  Physically,  the  scaling  factor  increases  the  initial  thrust  angle,  as 
seen  in  Figure  5.5  since  the  initial  thrust  angle  (j)  is  measured  up  from  the  horizontal  At,(0) 
axis.  In  the  case  with  circular  velocity  boundary  conditions,  the  initial  velocity  must  be 
greater  than  the  final  velocity  on  a  larger  circular  orbit.  By  raising  the  initial  thrust  angle, 
the  scaling  factor  decreases  the  initial  acceleration  in  the  tangential  direction,  reducing  the 
amount  of  velocity  to  be  removed  by  the  final  time.  In  the  case  with  =  1  instead  of 
=  0,  a  larger  initial  thrust  angle  is  needed  to  counter  the  gravity  while  raising  the  orbital 
altitude. 

If  the  two  velocity  components,  Equations  (5.18)  and  (5.19),  are  used  to  eliminate 
their  common  logarithmic  term,  the  following  relationship  is  obtained: 

-  x(o)^  ^  ^  ^m)2 m-j  ^  ^  ,  0  (531) 

The  initial  velocity  components  are  i(0)  =  0  and  y(0)  =  1  on  the  starting  circular  orbit. 
The  final  Cartesian  velocity  components  are  not  known  individually  at  t/,  and  they  may 
be  set  equal  to  zero  as  an  approximation  which  results  in  the  simplest  form  of  the  solution. 
This  has  been  found  to  be  the  most  useful  approximation,  but  others  have  been  tried  such 
as  x{tj)  =  0  and  y{tf)  =  1.  With  both  final  velocity  components  set  equal  to  zero,  the 
scaling  factor  may  be  introduced  to  produce  an  equation  in  a  single  variable  with  two 
parameters,  R  and  A,  Then,  Equation  (5.31)  becomes: 

-  ^  +  ^/(9  -  52  +  (9252  -  2^52)'  -  =  0  (5.32) 

2T. 
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The  scaling  factor  q  may  be  expressed  as  an  infinite  series  in  the  quantity  (l/A). 
After  solving  for  the  coefficients,  the  first  three  terms  in  this  series  solution  are: 
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1  R 
2A  4A2 


(5.33) 


Numerically,  it  has  been  found  that  truncation  after  the  third  term  is  adequate  in  this 
approximation  because  there  are  greater  sources  of  error  from  the  other  assumptions,  such 
as  zero  gravity.  Using  this  result,  the  initial  costates  are  given  by: 

‘  0-^+4^)'' 


Since  the  expansion  for  q  involves  the  quantity  (l/A),  it  should  not  be  used  if 
A  <  1,  If  A  <  1,  the  approximate  initial  costates  may  be  obtained  from  Equations  (5.26) 
and  (5.27),  which  do  not  include  the  scaling  factor  q.  The  time  of  flight,  is  not  as  sen¬ 
sitive  to  the  presence  of  the  gravity  term,  so  the  formula  for  tf  does  not  require  additional 
scaling. 

To  summarize.  Equations  (5.26)  and  (5.27)  give  the  approximate  initial  costate  val¬ 
ues,  and  Equation  (5.28)  is  the  approximate  time  of  flight.  If  A  >  1,  the  scaled  Equa¬ 
tions  (5.34)  and  (5.35)  should  be  used  to  improve  the  approximate  initial  costate  values. 
These  relationships  were  derived  using  the  assumptions  of  zero  gravity,  zero  mass-flow  rate, 
and  zero  final  velocity  components.  Even  using  these  assumptions,  the  results  lead  to  a 
good  initial  guess  for  the  associated  boundary  value  problem  over  the  parabolic  region  of 
the  costate  locus,  where  the  parameter  S  is  less  than  or  equal  to  one.  The  convergence 
sensitivity  of  this  model  will  be  discussed  at  the  end  of  the  chapter,  along  with  the  elliptic 
and  spiral  region  models. 
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Figure  5.6  Final  R  vs.  Orbital  Revolutions 


5,2  Elliptic  and  Spiral  Regions 

If  the  transfer  occurs  in  the  neighborhood  of  one-quarter  to  one-half  revolution,  the 
initial  costates  do  not  lie  near  the  parabolic  arc  or  the  point  (1,0).  In  this  case,  the  locus 
has  a  nearly  elliptical  shape  as  seen  in  Figure  5.4,  and  the  range  of  S  is  roughly  1  <  5  <  10. 
If  ^  =  2,  for  example,  5=1  corresponds  to  a  transfer  over  about  0.2  of  a  revolution,  and 
5  =  10  corresponds  to  a  transfer  over  about  3  revolutions.  As  R  increases,  the  range 
of  orbital  revolutions  narrows  between  lines  of  constant  5=1  and  5  =  10,  as  shown  in 
Figure  5.6.  The  well  known  Earth-Mars  transfer  example  given  in  Bryson  and  Ho  [7]  with 
R  =  1.525  lies  in  this  region  with  5  =  1.933,  and  has  roughly  0.3  revolutions.  This  example 
will  be  presented  in  Chapter  6. 

It  is  difficult  to  provide  a  physical  explanation  for  the  elliptical  form  of  the  locus  in 
this  region,  but  the  appearance  naturally  leads  one  to  try  the  equation  of  an  ellipse  as  a 
model.  In  the  parabolic  region,  the  thrust  term  in  the  equations  of  motion  is  dominant  over 
the  gravity  term.  In  the  spiral  region,  the  gravity  dominates  the  thrust.  Thus,  the  gravity 
or  thrust  may  be  taken  to  be  zero  as  a  limiting  case,  allowing  for  closed-form  solutions 
of  the  equations  for  the  states  and  costates.  The  zero-gravity  development  is  presented 
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Elliptic  Region  Angle 
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Figure  5.7  Elliptic  Region  Angles  for  R  =  1.1,5,100 

in  the  previous  section,  and  the  zero-thrust  development  is  in  the  literature  [26],  with  the 
results  described  later  in  this  section.  However,  the  thrust  and  gravity  terms  are  typically 
within  an  order  of  magnitude  of  each  other  in  the  elliptic  region.  Therefore,  neither  term 
may  be  considered  small  relative  to  the  other.  Several  attempts  were  made  to  develop  a 
perturbation  solution  for  this  region,  but  because  both  the  gravity  and  thrust  terms  must 
be  retained,  the  attempts  were  unsuccessful.  Since  there  is  no  closed  form  solution  of  the 
full  equations  for  the  states  and  costates,  we  approximate  the  initial  costates  in  the  elliptic 
region  with  the  equations  of  ellipses  through  curve  fitting  of  the  numerical  results,  using 
the  parameters  R  and  A. 

An  empirical  relationship  between  the  polar  angle  and  the  fourth  root  of  thrust 
acceleration  A  may  be  obtained  from  the  data  shown  in  Figure  5.7.  This  relationship  is 
used  in  a  linear  approximation  for  the  polar  angle  of  the  elliptic  region  to  construct  the 
complete  parametric  fit.  The  quarter  power  relationship  is  a  result  of  simple  trial  and  error 
to  produce  a  fairly  linear  graph.  Because  of  the  form  of  the  elliptic  model,  it  is  necessary  to 
relate  the  parameter  A  to  an  angle  along  the  elliptic  curve.  An  example  of  the  parametric 
fit  is  shown  in  Figure  5.8  for  R  —  1.525,  the  distance  to  Mars,  and  for  R  =  100.  The  model 
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Elliptic  Region  Fit 
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Figure  5.8  Elliptic  Fit  for  R  =  1.525,100 


accuracy  decreases  with  decreasing  R  and  A,  but  is  still  within  the  numerically  determined 
radius  of  convergence  using  the  modified  quasi-Newton  method  described  in  Chapter  4. 
The  elliptic  models  are  as  follows: 


A„(0) 


_ _  4(0-25)  ,  _ _ (r  ogN 

(0.498719 -0.811477i?(o  25))  ^  (0.279574  -  0.050554.R(-o  25))^  •  > 

-0.00091i?+ 0.4114 

1.0  -f  (0.00264i?  -f  0.616)  cos  +  0.0008124iZ  -|-  0.11876)  ^  ‘  ^ 

sin  -1-  0.03588  (5.38) 


At,(0)  =  0.165989  -  feo  cos 


(5.39) 


The  spiral  point  on  the  costate  locus  corresponds  to  problems  with  more  than  10 
revolutions,  and  it  is  well  known  [26]  that  the  optimal  thrust  direction  is  nearly  tangential 
to  the  path  through  most  of  the  transfer.  For  such  cases,  a  good  starting  guess  is  Au(0)  =  0, 
and  A„(0)  =  1.  The  time  of  flight  for  a  many-revolution  transfer  is  approximately  given 
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by  [26]: 


(5.40) 


To  solve  the  boundary  value  problem  using  the  previous  results,  the  parameter  S  is 
first  calculated  from  Equation  (5.25).  If  1,  then  Equations  (5.26),  (5.27)  and  (5.28) 
should  be  used  to  initialize  the  problem.  If  1  <  5  <  12,  then  Equations  (5.28)  and  (5.38)- 
(5.39)  provide  a  reasonable  set  of  initial  values  for  the  problem.  The  parameter  value  of 
5  =  12  is  used  to  extend  the  elliptic  costate  model  further  into  the  spiral  region.  If  S  >  12, 
the  initial  values  of  the  costates  will  be  close  to  Ai,(0)  =  0,  and  Av(0)  =  1.  The  approximate 
time  of  flight  for  this  case  is  given  by  Equation  (5.40). 

5.3  Convergence  Sensitivity 

To  measure  the  convergence  sensitivity  of  the  shooting  method  for  various  values  of 
S  and  a  plot  of  modified  quasi-Newton  method  iterations  is  presented  in  Figure  5.9. 
The  number  of  iterations  for  a  given  R  and  A  is  a  good  measure  of  the  success  of  the 
approximate  models,  since  closer  initial  values  require  fewer  iterations  for  convergence  to  a 
desired  error  tolerance  (as  defined  in  Appendix  A).  Generally,  the  sensitivity  of  the  system 
increases  with  decreasing  A  and  thus  increasing  5',  because  the  flight  times  become  much 
longer  for  a  given  R  value.  In  other  words,  the  “shots”  become  much  longer  in  the  shooting 
method.  The  initial  values  of  the  Lagrange  costates  tend  to  stay  within  the  neighborhood 
of  unity,  as  seen  in  Figure  5.4.  However,  the  flight  time  may  become  large  for  small  A, 
so  errors  in  the  time  models  tend  to  be  magnified  for  large  values  of  S  corresponding  to 
multiple  revolution  transfers.  This  explains  the  increased  system  sensitivity  for  large  5', 
but  the  modified  quasi-Newton  method  wiU  still  converge  using  the  models  for  the  time 
and  initial  costates  over  the  specified  ranges  of  A  and  R.  In  Figure  5.9,  the  S  parameter 
is  plotted  on  the  horizontal  axis  using  a  logarithmic  scale,  and  the  three  decades  shown 
correspond  to  the  three  regions  of  the  costate  loci  in  Figure  5.4.  There  are  peaks  in  the 
sensitivity  near  5*  =  1  and  S  =  12,  which  occur  in  the  transition  regions  between  the 
parabolic,  elliptic,  and  spiral  point  models  for  the  approximate  initial  costates.  Some  of 
the  curves  do  not  span  the  entire  range  of  5,  because  of  the  relationship  i?  =  +  1. 


5-17 


Convergence  Sensitivity  to  Initial  Costate  Model 


Figure  5.9  Convergence  Sensitivity  to  Initial  Costate  Model 


For  instance,  on  the  curve  of  constant  A  =  0.5,  Ji  =  201  when  S  =  20.  Since  this  value 
of  Jl  is  larger  than  the  range  of  the  models,  the  data  is  stopped  when  Ji  reaches  a  value 
of  100.  The  value  of  =  100  is  used  as  an  upper  limit  because  of  physical  considerations. 
In  an  Earth- centered  system,  the  lunar  orbit  is  at  roughly  60  DU,  based  on  the  radius  of 
the  Earth.  In  a  Sun-centered  system,  the  planet  Pluto  is  at  roughly  40  DU,  based  on  the 
radius  of  the  Earth’s  orbit  about  the  Sun.  Thus,  R  =  100  DU  is  more  than  adequate  to 
cover  the  range  of  currently  practical  missions. 

It  is  important  to  recognize  that  the  iteration  numbers  shown  in  Figure  5.9  are  based 
on  the  approximate  costate  models.  If  the  boundary  value  problems  represented  by  each 
point  on  the  curves  had  been  solved  by  slowly  varying  the  thrust  as  described  at  the 
beginning  of  this  chapter,  the  step  size  could  be  kept  small  enough  such  that  each  case 
could  be  done  in  about  5  iterations.  However,  that  technique  requires  a  slow  approach  from 
some  elusive  known  case  to  the  case  of  interest.  Using  the  approximate  costate  models 
provided  here  along  with  the  modified  quasi-Newton  method,  any  case  of  interest  may  be 
solved  directly  with  no  prior  knowledge.  Thus,  one  may  proceed  immediately  with  any 
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values  of  R  and  A  to  Figure  5.9  and  obtain  a  rough  idea  of  how  many  iterations  to  expect 
using  the  modified  quasi-Newton  method.  Without  the  initial  costate  approximations,  it 
would  be  unreasonable  to  expect  convergence  at  all  for  arbitrary  values  of  R  and  A,  unless 
one  is  very  fortunate. 

The  extreme  sensitivity  of  the  minimum-time  continuous  thrust  problem  to  the  initial 
costate  values  is  well  documented  [7,  10].  As  mentioned  earlier  regarding  the  continuation 
method  [12],  the  step  size  in  the  problem  parameter  A  must  be  constantly  decreased  to 
maintain  convergence  in  a  fixed  number  of  iterations  as  the  locus  moves  toward  the  center 
of  the  spiral  region,  where  A  approaches  zero.  This  decreasing  step  size  phenomenon  is  a 
direct  measurement  of  the  convergence  sensitivity,  but  this  also  depends  on  the  robustness 
of  the  search  technique.  Because  of  the  interdependence  of  the  model  accuracy  and  the 
capabilities  of  search  method,  both  the  model  and  the  modified  quasi-Newton  method  are 
used  to  produce  the  iteration  data  shown  in  Figure  5.9  instead  of  the  continuation  method. 
Clearly,  a  different  initial  costate  model  and  search  method  wiU  produce  different  iteration 
data.  Thus,  the  convergence  sensitivity  shown  in  Figure  5.9  is  only  intended  to  represent 
the  behavior  of  the  models  and  techniques  presented  in  this  dissertation. 

5-4  Optimal  Initial  Costate  Locus  under  the  KS  Transformation 

Figure  5.10  shows  the  optimal  initial  costate  locus  under  the  KS  transformation  for 
ii  =  2  with  m  =  0.  The  behavior  of  the  locus  is  similar  to  the  Cartesian  case  in  that 
the  spiral  region  corresponds  to  small  A,  many  revolution  transfers.  Also,  the  “top”  of 
the  curve  extends  upward  with  increasing  R.  The  entire  initial  costate  locus  is  stretched 
horizontally  compared  to  the  Cartesian  case.  The  focus  of  the  spiral  is  at  the  point  (2,0), 
which  may  be  found  by  letting  A  =  0  and  equating  the  KS  Hamiltonian  to  the  polar 
Hamiltonian  with  ficticious  time.  As  a  result,  the  initial  costate  should  be  taken  as 

2A<;(0)  for  the  purpose  of  initial  approximation.  Because  of  the  similarity  to  Figure  5.1,  one 
might  be  tempted  to  conclude  that  the  approximate  initial  costate  models  may  be  used  to 
initialize  the  problem  under  the  KS  transformation  in  all  cases.  However,  the  horizontal 
stretching  phenomenon  is  not  linear,  since  the  first  horizontal  axis  crossing  in  Figure  5.1 
is  not  exactly  twice  that  of  Figure  5.10,  even  though  the  spiral  point  is  located  at  (2,0). 
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KS  Costate  Locus 


Figure  5.10  Optimal  Initial  Costate  Locus  under  the  KS  Transformation 

Thus,  the  approximate  initial  costate  models  provide  only  a  rough  estimate  to  initialize  the 
problem  under  the  KS  transformation.  Convergence  may  still  be  achieved  in  some  cases, 
and  an  example  of  this  is  given  in  Chapter  6. 

5.5  Summary 

This  chapter  presents  an  original  approach  to  solving  the  minimum-time  orbital 
transfer  problem  under  continuous  thrust  by  modeling  the  initial  values  of  the  Lagrange 
costates.  To  accomplish  this,  numerical  results  from  many  cases  are  used  to  depict  graphi¬ 
cally  the  functional  relationship  between  the  initial  costates  and  the  parameters  R  and  A, 
Then,  analytical  and  empirical  means  are  used  to  provide  approximate  expressions  for  the 
initial  costates  and  time  of  flight.  The  success  of  the  approximate  models  is  examined  by 
presenting  the  number  of  modified  quasi- Newton  steps  required  for  convergence  over  the 
complete  range  of  R  and  A.  Again,  convergence  would  not  be  expected  at  all  for  arbitrary 
parameter  values  without  the  approximate  models,  since  there  is  no  other  information 
available  for  initialization.  With  the  approximate  models,  one  may  quickly  produce  a 
minimum-time  orbital  transfer  example  based  on  any  practical  spacecraft  design  and  mis¬ 
sion  requirement.  Finally,  a  locus  of  the  initial  costates  under  the  KS  transformation  is 
included  for  comparison. 
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VL  Numerical  Examples 


This  chapter  presents  a  series  of  numerical  examples  which  correspond  to  each  of  the 
earlier  theoretical  developments.  The  boundary  value  problems  are  solved  using  the  shoot¬ 
ing  method,  and  the  quasi-Newton  step  is  modified  with  the  dynamic  scaling  technique 
described  in  Chapter  4.  The  trajectories  are  propagated  with  an  integration  tolerance  of 
10“®,  and  an  error  tolerance  in  the  final  state  values  of  10“^  in  canonical  units.  A  discussion 
of  convergence  criteria  is  presented  in  Appendix  A. 

The  approximate  initial  costate  models  developed  in  Chapter  5  are  for  the  circle-to- 
circle,  coplanar  problem.  Therefore,  we  first  present  two  examples  of  coplanar,  circle-to- 
circle  problems,  to  demonstrate  the  parabolic  and  elliptic  models.  The  first  example  is  a 
military  application  with  very  low  thrust  [10],  and  the  second  is  the  well-known  Earth-to- 
Mars  example  given  by  Bryson  and  Ho  [7].  The  reader  should  note  that  the  approximate 
initial  costate  models  developed  in  Chapter  5  lead  to  convergence  for  all  values  of  A  and 
R  in  the  domain  of  intent  (0<i4<oo,l<i2<  100).  The  KS  transformation  is  used 
on  the  Bryson  and  Ho  example  for  comparison,  and  is  also  used  on  a  many-revolution 
circle-to-circle  problem,  which  uses  the  spiral  point  model.  This  last  example  may  be  used 
as  a  sample  truth  model  for  a  perturbation  technique  in  which  A  is  treated  as  a  small 
perturbation. 

The  approximate  initial  costate  models  may  also  by  used  on  problems  with  noncircu¬ 
lar  end  conditions  that  are  elliptical  or  hyperbolic.  Thus,  an  example  of  a  coplanar  transfer 
to  a  non-circular  end  condition  is  presented.  The  final  orbit  is  hyperbolic,  but  the  initial 
costates  stiU  provide  reasonable  starting  values,  and  convergence  is  achieved.  An  elliptical 
end  condition  case  is  presented  later  with  the  three-dimensional  examples. 

Continuing  away  from  the  coplanar,  circle-to-circle  case,  the  next  example  shows 
a  trajectory  from  one  circular  orbit  to  another  with  a  change  in  inclination.  Thus,  the 
trajectory  occupies  three  dimensions.  Convergence  is  stiU  achieved  using  the  initial  costate 
models.  An  example  with  the  largest  departure  from  the  coplanar,  circle-to-circle  case  is 
then  presented  with  an  inclination  change  combined  with  non-circular  end  conditions.  The 
final  orbit  in  this  example  is  both  polar  and  eUiptical. 
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Two  additional  examples  are  given  that  result  from  the  solution  of  multiple  boundary 
value  problems  with  varying  parameters.  In  the  first  example,  after  completing  many 
coplanar,  circle-to- circle  cases,  curves  of  time  per  revolution  are  plotted  for  a  large  range 
of  A,  with  different  values  of  R.  In  this  way,  it  is  possible  to  observe  the  behavior  of  the 
transfer  angle  as  a  function  of  R  and  A.  In  the  final  example,  a  plot  is  provided  that 
shows  the  existence  of  a  minimum-time  ascending  node  value  associated  with  i2,  A,  and 
inclination  in  the  noncoplanar,  circle- to- circle  case. 

6.1  Two-Dimensional  Problems 

6.1.1  Circle-to- Circle.  The  first  example  involves  a  spacecraft  design  that  has 
been  investigated  by  past  researchers  [10].  A  2400  kg  spacecraft  is  at  geosynchronous 
altitude  with  a  thrust  of  1.3  N  and  a  mass-flow  rate  of  —0.000069  kg/s.  In  canonical  units, 
the  problem  parameters  are  A  =  0.002425  and  rh  =  —0.000395.  Suppose  it  is  desired  to 
make  a  smaD  increase  in  altitude  in  minimum  time  to  avoid  a  ballistic  anti-satellite  weapon. 
If  J?  =  1.000336,  the  new  orbit  radius  will  increase  by  roughly  15  km.  In  this  case,  the 
quantity  S  =  ^  1)M  =  0.3722,  which  is  less  than  unity.  Thus,  the  optimal  initial 

costates  should  lie  near  the  parabolic  region  of  the  costate  locus,  and  the  approximate 
initial  values  are  given  by  Equations  (5.26)  and  (5.27).  The  approximate  time  of  flight  is 
given  by  Equation  (5.28).  Using  the  above  problem  parameters,  the  results  in  canonical 
units  are  given  in  Table  6.1. 


Table  6.1  AS  AT  Avoidance  Example 


iteration 

tf 

A„(0) 

A„(0) 

1 

0.7441893 

0.3720947 

0.1384544 

4 

0.7366198 

0.3395791 

0.1201369 

The  first  line  in  Table  6.1  is  the  set  of  initial  approximations  for  the  AS  AT  avoidance 
example.  The  first  iteration  is  considered  to  be  the  evaluation  of  the  approximate  models 
for  the  initial  costates  and  time  of  flight.  The  second  line  shows  the  final,  converged  values 
after  4  iterations  of  the  shooting  method.  The  converged  initial  costate  values  are  plotted 
in  Figure  5.4  as  a  point  on  the  costate  locus  which  is  identified  with  the  label  “ASAT.” 
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Figure  6.1  Iterative  Search  History  for  AS  AT  Avoidance  Example 


The  initial  costate  approximations  are  obtained  from  Equations  (5.26)  and  (5.27),  and  the 
approximate  time  of  flight  is  given  by  Equation  (5.28).  The  total  thrusting  time  is  roughly 
2.8  hours  in  physical  units.  In  this  case,  the  initial  approximations  for  the  costates  and 
flight  time  are  very  close  to  the  converged  values,  allowing  for  quick  convergence.  Figure  6.1 
shows  the  iterative  search  history  for  the  flight  time  and  initial  costates  from  the  initial 
guess  to  the  converged  values. 

Figure  6.2  compares  the  exact  control  angle  history  from  the  converged  case  with 
the  approximate  history  generated  by  the  initial  estimates.  Both  curves  pass  through  the 
zero  angle,  showing  the  switch  in  thrust  direction.  The  diflFerences  between  the  initial 
estimates  and  the  converged  values  are  most  evident  here,  since  the  costate  histories  are 
very  sensitive  to  the  initial  conditions. 

Figure  6.3  shows  a  comparison  of  flight  path  trajectories  using  the  control  law  from 
the  initial  estimates  and  the  converged  values.  In  this  case,  the  paths  are  almost  identical 
in  spite  of  the  differences  in  the  control  angle  histories.  Thus,  approximate  initial  costates 
obtained  assuming  zero  gravity,  zero  mass  flow  rate  and  zero  final  velocity  are  suflhcient  to 
achieve  an  almost  optimal  trajectory. 
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Figure  6.2  Control  Angle  History  for  ASAT  Avoidance  Example 


Figure  6.3  Optimal  Trajectory  for  ASAT  Avoidance  Example 


The  next  example  will  be  the  well  known  Earth-to-Mars  orbital  transfer  case  given 
in  Bryson  and  Ho  [7],  page  66,  where  R  =  1.525,  A  =  0.1405,  and  m  =  —0.07488  mass 
units  per  TU.  The  initial  costate  values  are  not  given  in  the  reference.  In  this  case,  the 
quantity  S  =  1.933,  so  1  <  5  <  10.  Thus,  the  optimal  initial  costates  should  lie  in 
the  elliptic  region  of  the  costate  locus,  and  the  approximate  initial  values  are  given  by 
Equations  (5.38)  and  (5.39).  Using  these  problem  parameters,  the  results  in  canonical 
units  are  given  in  Table  6.2. 


Table  6.2  Bryson  and  Ho  Example 


iteration 

A„(0) 

A.(0) 

1 

3.8660858 

0.4221273 

1.1093899 

7 

3.3192600 

0.4949228 

1.0785465 

The  first  line  in  Table  6.2  is  the  set  of  initial  approximations  for  the  Bryson  and 
Ho  example.  The  second  line  shows  the  final,  converged  values  after  7  iterations  of  the 
shooting  method.  The  converged  initial  costate  values  are  plotted  in  Figure  5.4  as  a  point 
on  the  costate  locus  which  is  identified  with  the  label  “Bryson-Ho.”  The  approximate  time 
of  flight  is  given  by  Equation  (5.28).  The  initial  approximations  for  the  costates  and  time 
of  flight  are  still  close  to  the  converged  values,  but  this  example  takes  more  iterations  for 
convergence  than  the  ASAT  example.  This  is  because  the  convergence  sensitivity  is  greater 
in  the  elliptic  region  than  in  the  parabolic  region,  and  the  initial  approximations  are  not  as 
close  as  in  the  ASAT  example.  However,  the  initial  approximations  are  aU  within  roughly 
15%  of  the  converged  values,  which  is  certainly  better  than  no  information  at  all. 

Figure  6.4  shows  the  iterative  search  history  for  the  flight  time  and  initial  costates 
in  the  Bryson  and  Ho  example  from  the  initial  guess  to  the  converged  values.  Figure  6.5 
compares  the  exact  control  angle  history  from  the  converged  case  with  the  approximate 
history  generated  by  the  first  and  third  iterations  of  the  initial  estimates.  All  curves  pass 
through  180  deg,  showing  the  switch  in  thrust  direction.  As  the  search  progresses,  the 
control  angle  history  approaches  the  optimal  control  angle  solution.  Again,  the  differences 
between  the  initial  estimates  and  the  converged  values  are  most  evident  here,  since  the 
costate  histories  are  very  sensitive  to  the  initial  conditions. 
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Figure  6.4  Iterative  Search  History  for  Bryson  and  Ho  Example 
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Figure  6.5  Control  Angle  History  for  Bryson  and  Ho  Example 
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Figure  6.6  Optimal  Trajectory  for  Bryson  and  Ho  Example 


Figure  6.6  shows  three  trajectories  corresponding  to  the  first,  third,  and  seventh 
iterations  in  the  shooting  method  to  solve  the  Bryson  and  Ho  example.  The  seventh 
iteration  is  considered  the  exact  solution  for  this  example,  as  described  in  Appendix  A. 
The  trajectory  from  the  first  iteration  overshoots  the  desired  final  radius  by  roughly  0.5  DU, 
and  the  third  is  much  closer.  The  “in  track”  error  along  the  spacecraft  velocity  direction  is 
dominant,  and  is  sensitive  to  changes  in  the  final  flight  time.  This  example  shows  typical 
behavior  in  which  the  final  flight  time  is  directly  related  to  the  final  radius,  assuming  the 
errors  are  small.  By  this,  we  mean  that  a  small  increase  in  final  time  over  that  of  the 
minimizing  path  wiU  result  in  a  small  increase  in  R.  Conversely,  a  small  decrease  from 
the  optimal  flight  time  will  result  in  a  small  decrease  in  i?.  It  is  the  observation  of  this 
direct  relationship  behavior  that  motivates  the  Jacobian  formulation  given  in  Chapter  4, 
Equation  (4.73),  in  which  the  final  time  is  used  instead  of  the  initial  value  of  A^. 


6.L2  KS  Transformation,  The  Bryson  and  Ho  example  will  now  be  examined 
under  the  KS  transformation,  for  comparison  with  the  previous  results  using  polar  coordi¬ 
nates  and  real  time.  Before  proceeding  with  the  full  example,  the  mass  will  be  held  constant 
to  allow  a  comparison  of  two  approaches,  with  and  without  a  separate  time  costate  A^. 
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Let  R  =  1.525,  A  =  0.1405,  and  m  =  0.  Given  these  parameter  values,  the  optimal 
fictitious  time  is  Sf  =  2.864381149.  If  is  not  used  in  the  constant-mass  case,  the  initial 
values  of  the  costates  are  as  follows:  (  Notice  A„,  =  2Au5,  and  A„j  is  from  .ff(O)  ) 

A„,(0)  =  -5.506657441 

A«,(0)  =  -2.957059105 

A,,(0)  =  -5.914118211 

A„,(0)  =  -12.948161442 

If  A(  is  used  in  the  constant-mass  case  and  Auj(O)  is  chosen  to  be  unity,  then  the 
initial  values  of  the  costates  are  as  follows,  and  Sf  is  unchanged: 

A„.(0)  =  1.000000000 

Au.(O)  =  0.536997105 

A,,(0)  =  1.073994210 

A„3(0)  =  2.351364974 

A,(0)  =  -0.181598367 

These  values  may  be  obtained  by  dividing  the  first  four  numbers  by  -5.506657441,  and 
also  At  =  (-1/5.506657441),  which  is  the  coefficient  of  the  Lagrangian  r  term.  Here,  the 
magnitude  of  the  largest  costate  has  been  reduced  by  roughly  a  factor  of  6,  making  the 
numerical  search  easier,  since  aU  of  the  numbers  are  closer  to  unity. 

Next,  the  complete  Bryson  and  Ho  example  is  presented  with  m  =  —0.07488.  The 
first  line  in  Table  6.3  is  the  set  of  initial  approximations.  These  initial  values  are  obtained 
from  the  approximate  models  of  Chapter  5,  and  are  the  same  values  that  are  used  in  the 
previous  presentation  of  this  example.  The  second  line  shows  the  final,  converged  values 
under  the  KS  transformation.  Figure  6.7  shows  the  iterative  search  history  for  the  flight 

Table  6.3  Bryson  and  Ho  Example  under  KS  Transformation 


iteration 

Au,(0) 

A«.(0) 

1 

3.8660858 

0.4221273 

2.2187798 

18 

2.7090520 

0.5371110 

2.3409681 

time  and  initial  costates  in  the  Bryson  and  Ho  example  under  the  KS  transformation. 
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Figure  6.7  Iterative  Search  History  for  Bryson  and  Ho  KS  Example 

from  the  initial  approximations  to  the  converged  values.  In  this  case,  convergence  took 
18  iterations  of  the  shooting  method.  The  number  of  iterations  is  larger  in  this  case  than 
in  the  previous  presentation  of  this  example,  mostly  driven  by  the  large  difference  in  the 
ficticious  time  and  the  real  time.  However,  the  approximate  initial  costate  and  flight  time 
models  were  not  derived  for  the  KS  transformation.  This  example  demonstrates  the  utility 
of  the  approximate  models  and  the  modified  quasi-Newton  method  in  solving  problems 
under  the  KS  transformation.  It  is  important  to  recognize  that  the  difference  between  t/ 
and  6f  will  increase  with  increasing  R,  so  convergence  may  not  be  achieved  for  large  values 
of  jR.  The  exact  control  angle  history  and  flight  path  are  the  same  as  shown  previously. 

As  a  final  example  using  the  KS  transformation,  the  initial  costate  values  and  ficti¬ 
cious  time  Sf  for  a  many-revolution  circle-to-cirde  transfer  are  presented  as  a  possible  test 
case  for  a  perturbation  technique.  If  the  thrust  is  taken  to  be  a  small  perturbation,  the 
resulting  minimum-time  trajectory  will  typically  be  a  many-revolution  spiral.  Thus,  the 
restdts  of  a  perturbation  analysis  may  be  checked  for  accuracy  against  this  example  for  the 
particular  values  of  R  and  A. 
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The  following  results  are  for  a  low  thrust  transfer  with  A  =  0.001,  R  — 2,  and  m  =  0. 
The  resulting  ficticious  time  is  Sj  =  215.5980372.  Again,  Xy^{0)  =  2Xu^(0)^  and  A^3(0)  is 
from  ^(0). 

A„^(0)  =  -945.557353841 

A«,(0)  =  17.418466660 

A^,(0)  =  34.836933320 

A^,(0)  =  -1999.696574002 

This  trajectory  takes  roughly  30  revolutions  to  complete,  and  resembles  a  simple 
spiral.  The  thrust  angle  never  exceeds  5  degrees  above  or  below  the  local  spacecraft  hori¬ 
zontal.  The  initial  values  of  the  costates  are  not  scaled  very  well,  as  the  largest  magnitude 
is  nearly  2000.  By  including  the  A^  costate,  aU  four  of  the  above  values  could  be  divided 
through  by  —945.557353841,  which  will  bring  them  all  nearer  to  the  interval  of  (—1,1). 
After  scaling  in  this  way,  A,,j(0)  =  -0.0368,  and  At;2(0)  =  2.1148,  which  is  close  the  spiral 
point  model  of  A,;^(0)  =  0,  and  A,,2(0)  =  2.  However,  for  simplicity,  it  may  not  be  desired 
to  include  A^  in  a  perturbation  technique.  The  reason  the  initial  costates  are  shown  with 
the  large  magnitudes  is  to  demonstrate  the  inherent  difficulty  brought  on  by  the  elimina¬ 
tion  of  At.  The  poorly  scaled  numbers  are  nowhere  near  the  interval  (  —  1, 1),  so  the  search 
becomes  more  dilficult  to  initialize.  If  it  is  required  to  eliminate  At  to  simplify  a  pertur¬ 
bation  technique,  the  poor  scaling  shown  above  must  be  retained  to  satisfy  ^^(0)  =  0,  as 
explained  in  Chapter  4. 

6.1.3  Circle-to- Hyperbola.  In  order  to  demonstrate  how  the  initial  costate  ap¬ 
proximations  may  be  used  for  a  problem  without  circular  end  conditions,  we  present  a  case 
in  which  the  final  conditions  are  hyperbolic.  The  approximations  may  also  be  used  if  the 
final  conditions  are  elliptical,  which  is  demonstrated  later  in  a  three-dimensional  exam¬ 
ple.  In  this  example,  a  spacecraft  starts  with  a  circular  orbit  of  unit  radius,  and  arrives 
at  i?  =  5  with  a  large  tangential  velocity.  The  final  orbit  is  hyperbolic  since  the  escape 
speed  is  0.6325  DU/TU,  and  the  spacecraft  has  a  much  larger  tangential  velocity  of  2.0 
DU/TU.  The  final  radial  velocity  is  zero,  so  the  spacecraft  arrives  at  the  periapsis  point 
of  the  hyperbola. 
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Table  6.4  HyperboHc  Example 


iteration  ; 

A„(0) 

A.(0) 

1 

10.7 

-0.039250 

1.226865 

10 

11.1 

-0.109450 

1.388047 

circ: 

8.1 

-0.130130 

1.388460 
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Figure  6.8  Iterative  Search  History  for  Hyperbolic  Example 


The  first  Hue  in  Table  6.4  is  the  set  of  initial  approximations  for  the  hyperbobc 
example.  The  second  line  shows  the  final,  converged  values  after  10  iterations  of  the 
shooting  method.  Table  6.4  also  shows  the  relationship  between  the  initial  costates  for 
the  circle-to-circle  case  and  the  hyperbolic  case.  The  third  fine  shows  what  the  converged 
initial  costates  and  flight  time  would  be  if  the  final  conditions  were  circular  at  =  5.  The 
biggest  difference  is  in  the  time  of  flight,  which  takes  longer  in  the  hyperbolic  case  to  reach 
the  final  desired  end  conditions. 

Figure  6.8  shows  the  iterative  search  history  for  the  flight  time  and  initial  costates 
in  the  hyperbobc  example  from  the  initial  guess  to  the  converged  values.  In  this  case, 
convergence  took  10  iterations  of  the  shooting  method. 

Figure  6.9  shows  the  control  angle  history  for  the  hyperbobc  example.  In  this  case, 
the  initial  thrust  angle  is  below  the  spacecraft  local  horizon.  Figure  6.10  shows  the  exact 
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trajectory  for  the  hyperbolic  example.  It  is  interesting  to  notice  that  the  optimal  solution 
takes  the  trajectory  beyond  R  =  5\n  order  to  give  the  spacecraft  “room”  in  terms  of  both 
space  and  time  to  accelerate  to  the  final  desired  velocity.  The  initial  and  final  thrusting 
times  are  indicated  on  the  trajectory,  to  show  that  the  final  hyperbolic  orbit  continues  on 
past  R  =  b. 

This  example  illustrates  the  idea  that  the  solution  trajectory  is  a  minimum- time  path 
to  a  complete  set  of  final  states,  not  just  a  maximum  radius.  Otherwise,  the  problem  would 
have  been  solved  as  the  spacecraft  first  reached  R  =  5.  However,  the  velocity  components 
were  not  yet  correct  at  that  time,  so  the  final  path  is  the  minimum-time  arc  to  the  final 
set  of  r,  u,  and  v. 

In  spite  of  the  highly  noncircular  final  conditions,  the  initial  costates  are  still  similar 
to  the  circle-to- circle  case,  as  shown  in  Table  6-4.  For  this  example,  convergence  is  achieved 
in  ten  iterations  of  the  shooting  method  using  the  approximate  initial  costate  model  based 
on  an  empirical  fit  to  the  circle-to-circle  case.  If  the  final  desired  end  conditions  are  for  an 
elliptical  path,  the  final  velocity  magnitude  would  lie  between  the  circular  and  hyperbolic 
cases.  Thus,  the  approximate  initial  costates  may  be  used  as  initial  guesses  for  elliptical 
end  conditions  as  well  as  for  the  circular  and  hyperbolic  end  conditions. 
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Figure  6.10  Optimal  Trajectory  for  Hyperbolic  Example 
6.2  Three-Dimensional  Problems 

To  show  the  relationships  between  the  two-  and  three-dimensional  cases,  a  two- 
dimensional  example  is  first  presented,  then  modified  to  include  an  inclination  change.  A 
third  example  is  presented  with  a  polar,  elliptical  end  condition.  This  example  uses  the 
converged  initial  costates  and  flight  time  from  the  second  example  as  initial  values.  All 
other  parameters  are  held  the  same. 

In  Figure  6.11,  an  optimal  transfer  is  shown  where  i?  =  5,  A  =  0.1,  and  m  =  —0.05. 
In  this  case,  the  parameter  S  takes  the  value  of  6.3246.  Accordingly,  the  initial  costates  lie 
in  the  elliptic  region  of  Figure  5.4.  To  solve  this  problem,  the  parameterized  elliptic  curve 
fits  are  used  to  provide  the  approximate  initial  costate  values. 

The  second  case  is  shown  in  Figure  6.12.  The  problem  parameters  are  the  same  as  in 
the  first  case,  but  the  desired  final  inclination  is  45  degrees,  and  the  desired  ascending  node 
is  zero  degrees  on  the  x  axis.  The  approximate  initial  values  for  and  A^  are  the  same 
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Figure  6.13  Flight  Path  for  Polar- Elliptical  Example 


as  in  the  two  dimensional  case,  and  A2(0)  and  A^y(O)  are  taken  to  be  zero.  Convergence  is 
achieved  in  14  iterations,  even  with  the  change  in  final  inclination. 

In  the  third  case,  the  goal  is  to  reach  a  polar,  elliptical  orbit  at  apoapse  where  the 
tangential  velocity  is  0.4  DU/TU,  which  is  less  than  the  circular  velocity  of  0.447  DU/TU. 
Convergence  is  achieved  in  18  iterations  using  the  three-dimensional  initial  values  from  the 
fourth  column  of  Table  6.4  with  no  further  modification.  Taking  large  steps  in  this  way  by 
using  continuation  combined  with  the  approximate  initial  costate  models  can  allow  for  the 
solution  of  a  wide  variety  of  orbit  transfer  problems  in  a  short  time.  In  this  case,  the  steps 
involve  inclination  changes  of  45  degrees  each  time.  The  converged  trajectory  is  shown  in 
Figure  6.13. 

Both  the  first  and  second  cases  are  solved  using  the  two-dimensional  approximations 
as  initial  values.  The  third  case  uses  the  results  of  the  second  case  as  initial  values.  The 
largest  difference  between  the  second  and  third  cases  is  in  the  flight  time.  Some  of  the  initial 
costates  in  the  third  case  are  actually  closer  to  the  approximate  two-dimensional  models 
than  those  in  the  second  case.  However,  the  difference  in  final  time  is  apparently  large 
enough  to  prevent  convergence  from  the  approximate  models  directly.  For  comparision, 
the  final,  converged  values  are  shown  in  Table  6.5  for  each  case. 

Figures  6.14,  6.15  and  6.16  show  the  iterative  search  histories  for  the  flight  time  and 
initial  costates  in  the  two-dimensional,  three-dimensional,  and  polar-eUiptical  examples 
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Table  6.5  Planar  vs.  Non-Coplanar  Example 


model 


2D 


W,i  =  45 


=  90 


A.(0) 

A.(0) 

A.(0) 

Ai(0) 

Ay(0) 

Ai(0) 

iterations  : 


1 

-0.1092950 

0 

-0.1092950 

1.1923589 

0 

12.6491106 


1 

-0.239259 

0 

-0.239259 

1.278794 

0 

10.045897 

12 


1 

-0.190287 

0.565053 

-0.190287 

1.281757 

0.443959 

10.985076 

14 


1 

-0.161830 

0.353082 

-0.161830 

1.550567 

0.827519 

14.437969 
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Figure  6.14  Iterative  Search  History  for  Two-Dimensional  Example 
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Figure  6.15  Iterative  Search  History  for  Three-Dimensional  Example 


3D  Polar-Bliptical  Example 
R,  =  5,A=0.1,i«90deg 


-^(0) 

-o-/^(0) 

-^(0) 


Figure  6.16  Iterative  Search  History  for  Polar- Elliptical  Example 
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from  the  initial  approximations  to  the  converged  values.  In  these  cases,  convergence  took 
12,  14,  and  18  iterations  of  the  shooting  method,  respectively. 

Figure  6.17  shows  that  the  behavior  of  the  thrust  angle  a  is  similar  between  the 
two-dimensional  and  three-dimensional  cases,  although  the  three-dimensional  case  has 
smoother  “corners,”  which  may  be  because  it  is  not  constrained  to  two  dimensions.  The 
thrust  angle  ^  starts  off  with  similar  behavior  between  the  three-dimensional  i  =  45°  case 
and  the  polar-eUiptical  example,  and  arrives  at  nearly  the  same  final  values  for  both  cases. 
Figure  6.18  shows  the  control  angle  histories  for  a  and  /?  in  the  polar-eUiptical  example. 
The  flattening  of  both  curves  at  the  end  of  the  trajectory  indicate  the  braking  maneuver 
needed  to  arrive  at  apoapse. 

In  the  two-dimensional  case,  the  costates  related  to  the  ^  direction  are  zero  through¬ 
out  the  transfer.  In  the  three-dimensional  case,  all  of  the  costates  vary  with  time.  The 
ascending  node  was  zero  for  this  example,  but  choosing  other  angles  will  result  in  diflFerent 
values  for  the  optimal  initial  costates  and  final  time.  The  best  agreement  between  the  two 
cases  is  for  Ay(0)  and  Ay(0),  which  may  be  approximated  with  Equations  (5.26)  and  (5.27), 
Equations  (5.34)  and  (5.35),  or  the  point  (0,1)  as  appropriate.  Convergence  is  achieved  in 
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Figure  6.18  Control  Angle  Histories  for  Polar- Elliptical  Example 


12  iterations  of  the  shooting  method  for  the  first  case,  14  iterations  for  the  second  case, 
and  18  iterations  in  the  third  case. 

6. 2 A  Time  per  Revolution.  Certain  trends  in  the  behavior  of  the  optimal  tra¬ 
jectories  become  evident  after  plotting  the  results  for  many  numerical  cases  [2].  As  the 
parameters  R  and  A  change,  the  time  of  flight  per  revolution  changes  as  well.  This  relation¬ 
ship  is  difficult  to  model  in  general  without  performing  numerical  integration.  However, 
the  maximum  value  may  be  modeled  as  a  function  of  R  with  a  fairly  simple  expression. 

Figure  6.19  shows  the  relationship  between  time  of  flight  and  transfer  angle  for  a 
range  of  values  of  A  and  R.  The  maximum  time  per  revolution  occurs  near  a  one-revolution 
transfer  for  each  value  of  R.  An  empirical  fit  for  the  one- revolution  case  is  as  follows: 

tj,  «  y/{R  -  0.6551)/0.01231  (6.1) 

A  plot  of  final  radius  against  the  time  of  flight  for  one  revolution  shows  a  roughly  parabolic 
shape,  as  seen  in  Figure  6.20. 

It  is  difficult  to  identify  a  physical  explanation  for  the  parabolic  relationship,  because 
the  path  is  a  powered  trajectory,  and  Keplerian  relationships  do  not  apply  directly  to 
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these  orbit  transfers.  The  empirical  relationship  in  Equation  (6.1)  is  found  by  assuming 
R  =  citj^  +  C2,  and  solving  for  the  unknown  coefficients  Ci  and  C2.  This  functional  form 
is  the  simplest  available  to  model  the  parabobc  shape  adequately.  A  numerical  example 
is  given  below  to  demonstrate  the  accuracy  and  utility  of  this  simple  model.  It  should  be 
noted  that  the  coefficients  in  Equation  (6.1)  are  found  using  a  small  number  of  data  points, 
so  this  is  at  best  a  rough  approximation.  The  data  points  are  found  manually  to  match 
the  requirement  of  exactly  one  revolution.  It  would  be  possible  to  modify  the  shooting 
method  to  include  the  unit  revolution  as  an  end  condition,  and  allow  the  thrust  magnitude 
to  vary  as  an  input.  However,  the  purpose  of  the  approximation  in  Equation  (6.1)  is  only 
to  provide  a  general  idea  of  the  relationship  between  the  time  of  flight  and  the  transfer 
angle  in  the  neighborhood  of  one  revolution. 

A  mission  planner  could  use  Equation  (6.1)  to  And  the  approximate  time  of  flight 
for  one  revolution  to  the  flnal  radius,  and  then  decide  whether  more  or  less  time  would 
meet  requirements.  An  estimate  of  the  required  thrust  may  be  obtained  by  using  this  time 
value  in  Equation  (5.28),  and  solving  for  the  parameter  A.  Similar  information  about  the 
relationship  between  time  of  flight  and  orbital  revolutions  is  presented  graphically  in  Ref¬ 
erence  [2].  As  a  numerical  example,  suppose  that  R  =  7  and  m  =  0.  Using  Equation  (6.1) 
and  Equation  (5.28),  we  have  tf^  «  22.7,  and  A  «  0.0466,  each  in  canonical  units.  The 
exact  values  are  =  22.9  and  A  =  0.0494  to  the  same  number  of  signiflcant  flgures.  Even 
as  a  rough  approximation,  Equation  (6.1)  gives  a  reasonable  value  of  the  optimal  flight 
time  for  a  one-revolution  transfer. 

6.2,2  Minimum  Time  vs.  Ascending  Node.  Figure  6.21  shows  the  variation  of 
time  of  flight  to  orbits  with  the  same  inclination  but  different  ascending  nodes.  As  the  flnal 
inclination  changes,  the  ascending  node  for  the  minimum  time  remains  relatively  constant. 
The  two  minima  of  each  locus  shown  in  Figure  6.21  are  of  equal  values,  and  they  represent 
trajectories  that  are  reflected  about  the  x,j/  plane.  The  values  of  the  ascending  node 
associated  with  these  minima  are  normally  spaced  at  tt  radians  apart.  This  is  because  the 
intersection  of  the  flnal  orbit  with  the  x^y  plane  is  a  line  segment,  and  the  two  minimum¬ 
time  paths  typically  arrive  at  locations  on  opposite  sides  from  the  origin  as  projected  onto 
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Figure  6.21  Time  of  Flight  to  Various  Ascending  Nodes 


the  line  segment.  It  is  difficult  to  model  the  minimum-time  node  as  a  function  of  R  and  A, 
due  to  extreme  sensitivity.  This  example  is  provided  to  show  a  qualitative  illustration  of 
a  typical  relationship  between  the  ascending  node  and  the  minimum-time  path  to  a  final 
orbit  with  a  given  inclination. 

6.3  Summary 

This  chapter  presents  a  series  of  numerical  examples  of  minimum-time  orbit  transfers 
under  continuous  thrust  for  two-  and  three-dimensional  end  conditions,  and  also  a  two- 
dimensional  problem  under  the  KS  transformation.  The  purpose  of  these  examples  is  to 
demonstrate  the  wide  range  of  applicability  of  the  approximate  initial  costate  models  and 
the  modified  quasi-Newton  method  developed  in  Chapter  5.  In  addition,  numerical  results 
from  many  converged  cases  are  presented  concerning  the  time  of  flight  for  one  revolution 
in  the  coplanar  circle-to-circle  case,  and  for  the  minimum-time  path  to  a  given  final  in¬ 
clination.  In  aU  cases,  the  approximate  intial  costate  relationships  result  in  convergence 
directly  or  by  continuation  with  a  large  step,  as  in  the  polar-elliptical  example. 
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VIL  Summary  and  Conclusions 


7. 1  Summary 

The  equations  of  motion  for  a  spacecraft  under  the  influence  of  gravity  and  continuous 
thrust  are  developed  using  several  coordinate  systems.  Euler-Lagrange  theory  is  then  used 
to  derive  the  optimal  control  law  and  diflFerential  equations  for  the  states  and  costates  in 
each  of  the  coordinate  systems.  Expressions  are  developed  for  the  approximate  optimal 
initial  costates  as  functions  of  R  and  A  for  the  minimum-time,  circle- to- circle  continuous- 
thrust  orbit  transfer  problem.  The  shooting  method  is  described  and  used  to  solve  the 
boundary  value  problem,  and  the  approximate  optimal  initial  costate  expressions  are  used 
for  the  initial  values.  A  dynamic  step  limiter  is  presented  which  improves  convergence 
in  the  shooting  method.  The  minimum-time  continuous-thrust  orbital  transfer  problem 
is  also  developed  under  the  Kustaanheimo-Stiefel  (KS)  transformation,  and  the  optimal 
initial  costates  are  presented  for  comparison.  Examples  are  provided  for  coplanar  and 
non-coplanar  orbital  transfers,  showing  the  utility  of  the  two-dimensional  approximations 
for  three-dimensional  problems  and  non-circular  end  conditions. 

7.2  Conclusions 

The  prime  motivation  for  this  work  is  that  there  are  no  models  available  in  the 
literature  that  provide  initial  costate  estimates  for  the  minimum-time,  continuous-thrust 
orbit  transfer  problem  as  functions  of  the  problem  parameters.  The  inherent  diflH.culty 
in  classical  optimization  methods  is  the  need  to  guess  the  initial  values  of  the  Lagrange 
multipliers.  The  expressions  and  techniques  developed  in  this  research  lead  to  convergence 
in  the  shooting  method  for  the  indicated  range  of  problem  parameters.  The  greatest 
advantage  to  this  approach  is  that  the  resulting  trajectories  satisfy  the  Euler-Lagrange 
equations,  so  optimality  is  guaranteed.  This  research  has  proven  successful,  in  that  it 
provides  a  reliable  means  to  determine  the  optimal  thrust  angle  history  from  the  shooting 
method  for  arbitrary  values  of  the  problem  parameters.  It  was  not  possible  to  do  this 
before  with  any  reliability,  because  the  only  other  option  is  to  rely  on  pure  guesswork  and 
good  fortune. 
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In  the  literature,  there  are  many  attempts  to  solve  approximations  of  the  optimal 
control  problem  by  making  simplifications  of  the  dynamics  equations  or  assumptions  about 
the  control  law.  However,  there  are  very  few  instances  where  researchers  provide  analysis 
for  the  problem  of  determining  the  initial  costates  [1,  13,  15,  18].  The  references  cited 
here  each  analyze  the  issue,  but  they  assume  impulsive  thrust  maneuvers  or  non-optimal 
constant  tangential  thrust.  There  is  no  information  from  these  sources  that  provides  any 
useful  information  to  solve  the  minimum-time,  continuous-thrust  problem  for  arbitrary 
values  of  the  parameters  R  and  A.  Thus,  the  research  presented  here  is  unique  in  that  it 
directly  addresses  the  optimal,  minimum-time  continuous-thrust  case  for  the  purpose  of 
modeling  the  optimal  initial  costates  based  on  the  problem  parameters. 

A  mission  planner  could  use  the  results  of  this  research  to  simplify  the  job  of  mission 
design  for  continuous-thrust  spacecraft.  Without  the  approximate  initial  costates,  it  is 
very  difficult  to  find  an  optimal  solution  to  match  desired  end  conditions.  If  low-thrust 
devices  are  used  on  military  spacecraft,  there  may  not  be  much  time  available  to  design  a 
tactical  orbital  maneuver.  This  adds  importance  to  achieving  convergence  to  the  optimal 
trajectory  as  quickly  as  possible. 

When  tracking  a  maneuvering  spacecraft,  it  is  usually  necessary  to  wait  for  the 
maneuver  to  be  completed  before  estimating  the  new  orbit.  Under  continuous  thrust, 
however,  the  maneuver  may  go  on  indefinitely.  In  order  to  track  such  a  spacecraft,  the 
estimation  process  must  include  some  sort  of  dynamics  model  that  takes  the  thrust  into 
account.  If  the  trajectory  of  the  thrusting  spacecraft  is  assumed  to  be  time  optimal,  the 
results  of  this  research  could  be  used  to  provide  a  reference  trajectory  for  the  estimation 
process. 

A  natural  extension  of  this  research  is  to  examine  the  optimal  initial  costates  for 
the  minimum-fuel  problem,  which  would  involve  throttling  and  coasting  arcs.  Also,  the 
behavior  of  the  optimal  initial  costates  might  be  modeled  for  continuous-thrust  Earth- 
Moon  trajectories,  where  the  Moon’s  gravity  is  included  in  the  equations  of  motion  [21]. 
The  variational  Hamiltonian  would  be  much  more  complex  in  the  restricted  three  body 
problem.  Another  extension  of  the  research  would  be  to  include  the  effects  of  uncertainty  in 
the  values  of  the  states,  and  to  examine  the  effect  on  the  costates  for  in-flight  computations. 


7-2 


There  are  many  possibilities  for  future  research  in  the  area  of  initial  costate  deter¬ 
mination  for  continuous-thrust  optimization.  Simplified  gravity  models  could  yield  closed 
form  solutions  for  the  equations  of  motion,  or  the  costate  equations.  Any  such  solution 
might  provide  an  analytical  approximation  for  the  initial  costates  when  the  accelerations 
due  to  gravity  and  thrust  are  of  nearly  the  same  magnitude.  Also,  the  KS  transforma¬ 
tion  development  and  examples  could  be  used  to  verify  pertubation  analysis  for  low-thrust 
problems. 
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Appendix  A,  Numerical  Solution  and  the  Shooting  Method 


To  solve  the  boundary  value  problems  presented  in  the  text,  the  differential  equations 
of  the  states  and  costates  must  be  numerically  integrated  from  the  initial  time  to  the  final 
time.  The  final  conditions  will  depend  on  the  choices  for  the  initial  conditions.  If  the 
final  conditions  are  incorrect,  then  the  initial  conditions  must  be  adjusted  according  to 
Equation  (4.73).  The  partial  derivatives  in  the  Jacobian  matrix  may  be  found  numerically, 
through  the  shooting  method.  First,  a  reference  trajectory  is  propagated  by  numerically 
integrating  the  differential  equations  from  the  initial  time  to  the  final  time.  Then,  small 
perturbations  are  made  in  each  of  the  unknown  initial  conditions,  and  the  trajectory  is 
propagated  again  individually  for  each  of  the  perturbations.  The  relative  errors  at  the 
end  time  are  collected  and  used  with  the  perturbation  magnitudes  to  calculate  one-sided 
approximations  of  the  partial  derivatives  in  the  Jacobian  matrix.  The  corrections  to  the 
initial  conditions  are  then  given  by  Equation  (4.73). 

Because  of  the  nonbnearities  in  the  differential  equations,  the  errors  in  the  final 
conditions  may  grow,  shrink  or  stay  the  same  as  the  shooting  process  is  repeated.  With 
good  initial  costate  models,  and  the  modified  Newton  step  described  in  Chapter  4,  the  final 
errors  wiU  normally  get  smaller  with  each  iteration  of  the  shooting  method.  Figure  A.l 
is  a  flowchart  for  the  shooting  method.  Many  sources  of  computer  code  are  available 
to  accomplish  the  tasks  shown  in  this  flowchart  [16].  The  code  is  available  in  the  form 
of  subroutines  in  various  computer  languages.  The  coding  for  this  research  was  done  in 
Borland  Turbo  Pascal  for  Windows,  version  7.0,  and  the  computations  were  performed  on 
a  Pentium-90  based  personal  computer. 

The  trajectories  in  aU  of  the  examples  in  this  dissertation  are  propagated  with  an 
integration  tolerance  of  10“®,  and  an  error  tolerance  in  the  final  state  values  of  10“^  in 
canonical  units.  This  defines  the  “exact”  solutions  referenced  in  the  text,  regardless  of 
which  iteration  in  the  shooting  method  results  in  the  required  error  tolerance.  Borland 
Turbo  Pascal  can  support  “extended”  precision,  which  gives  a  floating  point  variable  an 
accuracy  of  19  digits  past  the  decimal  point.  This  is  much  more  precision  than  10“®, 
but  it  is  not  necessary  to  exercise  the  full  capability  of  the  machine  and  programming 
language.  Based  on  numerical  experience,  convergence  always  occurs,  once  the  final  errors 
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STEP  1: 
READ  INPUT 


Figure  A.l  Flowchart  for  the  Shooting  Method 


are  reduced  to  the  order  of  10"^  or  10"^  in  canonical  units.  The  error  tolerance  of  10”'^ 
is  three  or  four  orders  of  magnitude  more  precise  than  necessary  to  indicate  convergence, 
and  thus  provides  confidence  that  the  minimizing  path  has  been  found  for  the  particular 
problem.  If  aU  19  digits  are  used,  the  integration  time  becomes  much  larger  per  trajectory. 
Since  the  shooting  method  requires  many  trajectories  to  be  calculated  during  an  iterative 
search,  the  increase  in  computing  time  is  multiplied,  with  no  great  benefit. 
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